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ABSTRACT 

Lewis,  Karl  Harcourt.   Ph.D.,  Purdue  University,  June  I967. 
"Analysis  of  Concrete  Slabs  on  Ground  and  Subjected  to  Warping  and 
Moving  Loads."  Major  Professor:  M.  E.  Harr. 

A  theory  has  been  developed  whereby  stresses  and  deflections  could 
be  calculated  for  a  series  of  rectangular  slabs  lying  on  a  viscoelastic 
foundation  and  subjected  to  a  moving  load.   The  stresses  and  deflections 
are  caused  by  the  weight  of  the  slab,  the  moving  concentrated  load,  and 
by  linear  temperature  (or  moisture)  variations  that  cause  sufficient 
warping  so  that  the  slab  is  only  partially  supported  by  its  foundation. 

Part  I  considers  the  problem  of  partial  support  caused  by  warping 
(in  this  thesis,  the  temperature  at  the  top  of  the  slab  is  smaller  than 
that  at  the  bottom),  while  Part  II  concerns  itself  with  the  effect  of  a 
reduction  in  subgrade  support  over  some  narrow  region. 

In  both  parts,  the  support  conditions  were  simulated  by  a  Kelvin 
viscoelastic  model,  and  zones  (which  depended  on  the  value  of  subgrade 
reaction)  were  set  up  so  that  the  solutions  to  the  governing  differential 
equations  could  be  reduced  to  a  set  of  simultaneous  algebraic  equations. 
For  the  problem  studied  in  Part  I,  the  resulting  simultaneous  equations 
were  non-linear  and  a  method  of  functional  iteration  equivalent  to  an 
N-diraensional  Newton's  Method  was  used.   In  the  case  of  Part  II,  the 
equations  were  linear  and  the  Method  of  Crout  Reduction  was  used.   In 
both  parts,  the  equations  were  solved  with  the  aid  of  an  IBM  7090 
digital  computer  using  a  Fortran  source  program. 
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It  was  found  that  when  partial  support  due  to  warping  exists,  the 
tensile  stress  in  the  slab  can  increase  with  increasing  velocity  of  load. 
Moreover,  the  maximum  deflection  (downward)  need  not  occur  when  the 
velocity  of  the  load  is  equal  to  zero.  The  reduction  in  subgrade  support 
over  a  narrow  region  (approximately  8  feet  or  less)  does  lead  to  deflec- 
tions and  stresses  which  are  higher  than  those  calculated  using  the 
initial  value  of  subgrade  reaction.   This  is  particularly  evident  when 
the  load  is  over  the  region  of  reduced  subgrade  reaction. 


INTRODUCTION 

Improvements  in  the  performance  of  concrete  pavements  for  highways 
and  airports  have  for  many  years  been  the  cause  of  concern  among  Civil 
Engineers  interested  in  these  problems.  This  concern  has  not  only  been 
due  to  the  ever  increasing  cost  of  construction  and  maintenance,  but 
also  to  the  preponderance  of  cracks  which  exist  in  pavements. 

The  development  of  these  cracks  depends  on  several  factors,  (such 
as  type  of  subgrade,  deterioration  of  the  concrete,  temperature  effects 
and  load)  (68)*  and,  although  they  may  or  may  not  represent  a  failed 
condition,  they  do  indicate  deficiencies  in  the  analysis,  design  or 
construction  of  highway  pavements.   Cracks  not  only  detract  from  the 
general  appearance  of  a  highway,  but  often  they  lead  to  driving  discomfort 
and  pavement  deterioration;  especially  when  load  transfer  and  subgrade 
support  are  lost. 

As  far  as  the  construction  of  concrete  pavements  is  concerned, 
advancement  in  the  development  of  the  principles  of  soil  mechanics  has 
made  great  strides  in  limiting  the  use  of  poor  subgrades  and  has 
virtually  eliminated  the  use  of  soils  susceptible  to  frost  action  and 
pumping.   In  addition,  specifications  (1,  63)  geared  at  promoting  the 
use  of  sound  concrete  and  the  practice  of  good  construction  methods,  have 
become  widespread. 


*  Figures  in  parenthesis  refer  to  references  in  Bibliography. 


Along  with  refinements  in  construction  procedures,  methods  of 
monitoring  the  factors  influencing  concrete  pavements  have  been  developed 
(21,  22),  and  strains  as  well  as  deflections  of  pavements  can  be  determined 
for  both  static  and  dynamic  loadings  (16,  2k,    5U).   Temperature  differences 
between  the  upper  and  lower  faces  of  the  slab  can  be  measured  with  the 
use  of  thermo-couples  (67),  and  variations  in  moisture  content  can  be 
established  from  the  dielectric  properties  of  concrete  (3). 

Today,  the  problem  of  material  weakness  has  been  reduced,  if  not 
solved,  and  significant  progress  has  been  made  in  the  ability  to  build 
pavements  and  to  measure  their  properties  which  seem  important;  however, 
only  little  has  been  done  to  upgrade  the  methods  of  analysis.   Current 
procedures  for  designing  and  evaluating  pavements  (68)  are  still  based 
on  static  loads  and  except  for  introducing  equivalent  static  loadings  (2), 
they  do  not  account  for  the  dynamic  response  of  the  pavement  to  moving 
loads.  Moreover,  the  pavement  is  usually  assumed  to  be  fully  supported 
at  all  times  even  though  this  often  is  not  the  case  (lb,  29,  30). 

The  present  trend  in  transportation  is  towards  heavier  loads  and 
higher  speeds  and  consequently  there  is  an  increasing  need  to  be  able 
to  predict  performance.   This  can  only  be  done  when  more  is  known  about 
the  input  and  causal  factors  of  traffic  loads,  curling,  and  the  loss  of 
subgrade  support.   It  is  to  help  satisfy  this  need,  that  this  thesis  is 
dedicated. 


REVIEW  OF  LITERATURE 

In  the  early  stages  of  development,  pavement  design  and  road 
construction  consisted  of  a  collection  of  rule  of  thumb  procedures  based 
on  empirical  knowledge.  As  early  as  1909,    test  roads  were  being  used  to 
determine  the  best  type  of  pavement  for  the  prevailing  traffic  (^9)- 
Between  1920  and  19^0,  much  work  was  done  on  the  classification  of  soils 
(5,  57)  and  Highway  Engineers  were  able  to  correlate  pavement  performance 
with  subgrade  types. 

With  the  advent  of  World  War  II,  engineers  were  faced  with  the 
problem  of  designing  pavements  for  greater  wheel  loads  than  previously 
thought  necessary.   They  realized  that  they  could  not  afford  the  long 
period  of  time  which  was  required  to  develop  design  procedures  based  on 
past  experiences  and  thus  sought  more  rational  basis  of  design.   This 
search  ultimately  led  to  procedures  which  form  the  core  of  pavement 
design  today. 

Static  Load  Solutions 
In  1884,  Hertz  (19)  first  published  a  solution  to  the  problem  of 
an  elastic  plate  on  a  Winkler-type  foundation  {66);   however,  it  was  not 
until  the  advent  of  Westergaard's  solution  in  1926  that  a  highway  pave- 
ment was  treated  in  this  manner  and  the  problem  was  approached  from  a 
purely  theoretical  point  of  view.   Today,  this  work  still  forms  the 
basis  of  the  analytical  bent  in  pavement  design. 


Westergaard  (6U)  solved  the  problem  of  a  slab  fully  supported  on  a 
Winkler  foundation  and  subjected  to  static  loads  applied  at  the  interior, 
free  edges  and  corners  of  the  slab.   Later,  Kelley  (32),  Spangler  (56), 
Pickett  (1+5)  and  Westergaard  (65)  himself  extended  these  original 
solutions  to  account  for  linear  temperature  variations.   Thomlinson  (58) 
in  I9U0  further  modified  Westergaard' s  solution  by  assuming  a  simple 
harmonic  temperature  variation  which  resulted  in  curved  temperature 
variations  throughout  the  depth  of  the  slab.   In  1957,  Freudenthal  and 
Lorsch  (15)  used  the  three  fundamental  models  (Maxwell,  Kelvin  and 
Standard  Solid)  to  study  the  problem  of  an  infinite  beam  on  a  visco- 
elastic  foundation;  and  in  I958,  Hoskin  and  Lee  (26)  solved  the  problem 
of  an  infinite  plate  on  a  viscoelastic  foundation. 

The  foregoing  solutions  neglect  the  shearing  forces  generated  at 
the  pavement-base  interface.   Several  mechanisms  have  been  offered  to 
account  for  this  effect.   Filonenko-Borodich  (13)  considered  a  set  of 
springs  held  together  by  a  membrane,  while  Schiel  (55)  took  a  fluid  which 
exhibited  surface  tension  as  his  soil  model.   Pasternak  (U3)  and  Kerr  (3*0 
on  the  other  hand,  considered  a  beam  of  unit  depth  resting  on  a  bed  of 
interrelated  springs  as  their  foundation,  and  Pister  and  Williams  (hk) 
used  the  shear  interactions  suggested  by  Reissner  (52).   Perhaps  the 
most  realistic  approach  offered  to  account  for  shear  in  an  elastic  base 
is  that  given  by  Klubin.   Klubin  (35)  expressed  the  pavement  reaction 
by  an  infinite  series  of  Tschebyschef f  polynomials  which  also  accounts 
for  the  two  elastic  constants  (Modulus  of  Elasticity  and  Poisson's  Ratio). 

It  should  be  noted  that  Klubin' s  solution  is  in  fact  an  elastic 
solution  while  the  others  which  impose  Winkler  assumptions  are  not.   In. 
spite  of  this,  since  the  Winkler  foundation  affords  a  much  simpler 


analysis  and  generally  gives  good  agreement  with  field  data,  it  will  no 
doubt  remain  popular. 

All  of  the  above  analyses  are  based  on  the  assumption  that  the  slab 
maintains  contact  with  its  support  at  all  points  and  at  all  times. 
Experimental  and  field  studies  (l6,  29,  30)  have  shown  this  assumption 
to  be  seriously  in  error  and  a  few  investigators  have  accounted  for  the 
effects  of  only  partial  support.   In  1959,  Leonards  and  Harr  (37)  solved 
the  problem  of  a  partially  supported  slab  on  a  Winkler  foundation  for 
linear  temperature  and/or  moisture  variations.   Later,  Reddy,  Leonards, 
and  Harr  (51)  extended  this  analysis  by  introducing  non-linear  temperature 
variations  as  well  as  a  viscoelastic  foundation.   In  all  of  these 
analytical  procedures,  only  the  case  of  symmetrical,  statically  applied 
loads  were  considered. 

Dynamic  Load  Solutions 
For  many  years,  the  problem  of  determining  the  stresses  and 
deflections  in  a  vibrating  plate  has  been  of  interest  primarily  to 
mathematicians.   Raleigh  (k8)    and  Lamb  (36)  using  the  classical  beam 
theory  developed  by  Euler  (ll)  and  Bernoulli  (k) ,    studied  several 
problems  dealing  with  the  vibration  of  bars,  membranes  and  plates. 
Later,  Ritz  (53)  elaborated  on  this  work  and  made  significant  contribu- 
tions toward  the  study  of  vibrating  rectangular  plates. 

Recently,  engineers  have  felt  the  need  to  account  for  the  dynamic 
response  of  pavements  and  several  solutions  have  evolved.   Pioneering 
these  solutions  was  the  work  of  Timoshenko  (60),  Hovey  (28)  and  Ludwig  (39) 
in  their  studies  of  the  dynamics  of  rails  subjected  to  moving  loads. 
In  19^3 »  Dorr  (10)  using  Fourier  integrals  extended  the  idea  to  a  beam, 


but  in  all  of  these  solutions  the  foundation  was  represented  by  a 
Winkler  model  which  exhibited  no  viscous  effects. 

Later  in  1953,  Kenney  (33)  added  the  effects  of  linear  damping  and 
by  means  of  the  method  of  undetermined  coefficients  was  able  to  examine 
the  relationship  between  deflections  and  critical  velocity.   Kenney 
found  that  at  a  constant  velocity  ratio  (velocity/critical  velocity) 
the  deflection  amplification  factor  (ratio  of  deflection  behind  moving 
load  to  that  under  static  load)  increased  as  the  damping  ratio  (damping 
coefficient/critical  damping  coefficient)  decreased.  For  heavy  damping 
(damping  ratio  =5r  1.5),  the  deflection  amplification  decreased  with 
increasing  velocity,  however  at  low  damping  (damping  ratio  t»    C.3)  this 
was  only  partially  true.   In  the  latter  case,  the  deflection  amplification 
increased  until  about  the  point  of  critical  velocity  (becoming  infinite 
for  a  damping  ratio  =0)  and  then  decreased.  In  the  case  of  the  moment 
amplification  factor  (ratio  of  maximum  moment  under  moving  load  to  that 
under  static  load),  a  similar  trend  was  experienced;  however,  at  low 
damping  ratios,  the  maximum  moment  occurred  at  velocities  slightly  below 
the  critical  value.  Mathews  (40)  in  1958  also  studied  a  similar  problem 
using  Fourier  transforms.   Crandall  (9),  recognizing  the  contributions 
of  Timoshenko  (60)  and  Mindlin  (41)  on  shear  deformation,  extended 
Kenney' s  solution  by  replacing  the  Bernoulli-Euler  beam  with  the  Timoshenko 
model.   Crandall  found  that  whereas  the  former  model  had  a  single  resonant 
frequency  and  a  single  critical  velocity,  the  latter  had  three  of  each. 

In  a  recent  work,  Thompson  (59)  elaborated  on  Kenney' s  solution  and 
showed  that  the  solutions  for  the  deflections  fall  intq  three  distinct 
regimes.  Specifically,  Thompson  showed  that  these  regimes  were-defined 
by  the  value  of  the  discriminant  of  the  fourth  order  characteristic 


equation.   If  the  value  of  the  discriminant  as  defined  by 

A  -  16  [  h   (1-e2)  98  -  (8  -  36c2  +  27e^)  9k  +  k\  (l) 

where   e  =  the  damping  ratio 

and    9  =  the  velocity  ratio 

is  greater  than  zero,  the  characteristic  equation  has  no  real  roots. 

If  A  is  equal  to  zero,  there  is  one,  real,  double  root;  and  if  A  is 

less  than  zero  there  are  two,  real,  unequal  roots. 

The  solutions  presented  by  Thompson  demonstrate  that  at  static 
conditions,  the  deflection  curve  is  symmetrical  (with  maximum  deflection 
occurring  under  the  load);  but  as  velocity  increases,  the  point  of 
maximum  deflection  falls  further  and  further  behind  the  load.   Upward 
deflections  occur  behind  the  moving  load  when  A  is  greater  than  zero, 
however  when  A  is  less  than  zero,  there  are  no  upward  deflections  (i.e. 
the  deflected  surface  does  not  intersect  the  axis  of  zero  deflection, 
but  simply  approaches  it  at  some  great  distance  behind  the  load). 

Several  investigators  have  also  considered  the  road  loading  system 
as  an  interaction  between  two  major  components  which  are  interdependent, 
namely  vehicles  and  roads.   Fabian,  Clark  and  Hutchinson  (12)  examined 
the  elements  of  each  component  and  developed  some  basic  mathematical 
models.  Analysis  of  their  vehicle  sub-system  showed  that  the  magnitude 
of  the  dynamic  load  is  a  function  of  vehicle  dynamic  properties  and 
apparent  road  profile,  and  may  in  fact  be  significantly  greater  than 
the  static  load.  Quinn  and  De  Vries  (h6)   used  an  experimental  procedure 
to  determine  the  highway  and  vehicle  characteristics  and  showed  that  if 
these  quantities  are  known,  the  reaction  of  a  vehicle  on  a  highway  can 
be  predicted. 
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Recognizing  the  difficulty  in  determining  appropriate  characteriza- 
tions for  both  highway  and  vehicle,  Sanborn  (5*0>  in  a  recent  work  used 
a  more  direct  approach.  He  simultaneously  measured  the  dynamic  load 
applied  to  a  pavement  by  a  moving  vehicle,  and  the  dynamic  stresses  and 
displacements  induced  by  that  load  in  the  pavement  and  its  subgrade. 
The  displacements  and  stresses  within  the  pavements  studied  varied  with 
velocity,  but  no  definite  pattern  was  established.  However,  dynamic 
variations  of  load,  for  the  load  and  range  of  velocity  covered  (1*500  lbs. 
and  0  to  60  mph  respectively),  appeared  to  have  little  effect  on  pave- 
ment response. 

In  general,  most  of  the  analytical  studies  reported  in  the 
literature  on  moving  loads  on  pavements  (68),  do  not  account  for  the 
interaction  between  vehicle  and  road.   For  the  most  part,  these  reports 
show  that  for  fully  supported  pavements,  the  lower  the  velocity,  the 
greater  the  deflections  and  stresses  In  the  slab.  Measurements  of  dis- 
placement (21,  23,  2^)  generally  support  this  finding;  however,  signif- 
icant data  (16)  does  exist  to  suggest  the  possibility  that  displacements 
may  increase  at  velocities  greater  than  30  to  kO  mph.  Also  to  be  noted 
is  the  fact  that  tests  on  tire  forces  on  pavements  {hi)   do  indicate 
that  the  "average"  force  produced  on  the  pavement  by  a  moving  load, 
increases  with  velocity. 

Obviously,  there  are  several  interdependent  factors  such  as  the 
vehicle's  suspension  system,  tire  pressures  and  road  profile  which  should 
be  included  in  the  analysis  of  pavements;  however  at  the  present  time, 
the  interaction  between  vehicle  and  pavement  is  not  well  understood.   In 
spite  of  this,  further  insight  into  the  highway  problem  can  be  gained  by 
examining  some  of  the  simplifying  assumptions  which  appear  in  current 


analyses  of  highway  pavements.   Thus,  it  is  the  object  of  this  thesis 
to  obtain  by  analytical  means  the  stresses  and  deflections  in  a  partially 
supported  concrete  slab  when  subjected  to  a  moving  load.  Such  partial 
support  may  be  caused  by  temperature  and  moisture  gradients  or  by  the 
weakening  and  partial  or  complete  loss  of  subgrade. 
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FORMULATION  OF  PROBLEM 

In  the  first  part  of  this  thesis,  the  problem  is  considered  of  a 
slab  lying  on  a  foundation  while  subjected  to  moving  loads  as  well  as 
to  temperature  (and/or  moisture)  gradients  which  would  cause  upward 
curling.*  The  second  part  will  investigate  conditions  where  a  weakening 
of  the  subgrade  has  developed.  The  conditions  in  Part  II  may  exist 
where  (a)  water  has  infiltrated  under  the  pavement  from  the  side  of  the 
roadway;  or  (b)  where  a  joint  has  opened  or  where  a  crack  has  occurred 
in  a  warped  pavement  to  such  an  extent  as  to  enable  it  to  regain 
complete  contact  with  the  foundation.  When  such  an  opening  exists,  water 
may  infiltrate  through  the  surface  of  the  pavement.   Under  both  conditions 
(a)  and  (b),  the  infiltration  may  weaken  the  subgrade  and  may  eventually 
lead  to  pumping.  The  following  assumptions  are  made  for  both  parts  in 
order  to  render  them  tractable: 


*  Upward  curling  is  understood  to  exist  when  temperature  gradients, 
cause  the  ends  of  the  slab  to  rise  vertically. 
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ASSUMPTIONS 

1.  A  highway  pavement  can  be  represented  by  an  array  of  rectangular 
plates. 

2.  The  usual  assumptions  of  plate  theory  hold,  namely:   the  plate  is 
homogeneous,  isotropic  and  obeys  Hooke's  law;  deflections  are  small 
in  comparison  to  thickness;  plane  cross  sections  normal  to  the 
middle  plane  in  the  unstressed  condition  remain  normal  to  this 
surface  after  bending;  the  effects  of  rotatory  inertia  and  shear 
deformation  can  be  neglected. 

3.  The  highway  base  material  acts  like  a  set  of  linear  viscoelastic 
elements.   The  inertia  of  the  material  is  neglected. 

k.     External  forces  acting  on  the  plate  are  those  due  to  a  constant- 
velocity  line  load  and  gravity. 

5-   The  plates  are  subjected  to  changes  in  temperature  (and/or  moisture) 
which  vary  linearly  with  depth.  The  variation  in  temperature  is 
constant  on  all  planes  parallel  to  the  upper  and  lower  plate  surfaces 
and  is  independent  of  time. 
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PART  I 


In  addition  to  the  foregoing  assumptions,  further  assumptions 
regarding  the  interaction  between  slabs  are  required.   For  the  case 
studied  here,  it  is  assumed  that  no  bending  moment  exists  between  slabs 
(this  will  be  discussed  later).   However,  some  load  transfer  can  be 
accounted  for  by  specifying  a  shearing  force  between  slabs,  equivalent 
to  that  existing  in  an  infinite  slab. 

In  general,  the  governing  differential  equation  describing  the 
pavement  section  illustrated  in  Figure  1  can  be  expressed  as  follows  (h2)  : 

tv  /Q  w   ,  ,_   d  w    ,  d  w   .  ,     d  w    /         »     ,        .N   ,  N 
D  (— jj-  +  2    2   2  +  — ^  )  +  PH  —  =  q(x   yv    t)  -  p(x  ,  y  ,  t)   (2) 

dx.      dx.  dy.    dy.         dt       l 
where  D  =  flexural  rigidity  of  the  slab 

w  =  mid-plane  deflection  of  the  slab  (positive  downward) 
x. ,  y.   =  fixed  coordinates 

p  =  density  of  the  slab 
H  =  slab  thickness 
q  =  surface  loading 
p  ■  foundation  reaction 
t  =»  time 

Using  Assumption  3>  the  foundation  reaction,  p,  may  be  expressed  as 

Sw 
(see  Figure  lb):   p(x.  ,  y. ,  t)  =>  C  r—  +  Kw,  where  C  is  the  damping 

coefficient  and  K  is  the  modulus  of  subgrade  reaction.  Applying  the 

loading  conditions  implied  by  Assumption  k,   and  assuming  that  the 
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Ik 


deflection  does  not  vary  in  the  y.  direction,  Equation  2,  for  a  constant 
cross  section  of  pavement,  becomes 


DS\+P„&  +  C|H  +  K„     +P(     , 

ox.       dt 

in  which  q  ■  unit  weight  of  slab 


(3) 


and       P  =  the  moving  line  load 

Employing  the  transformation,  x  ■  (x.  -  vt),  the  equation  is  seen 
to  reduce  to  a  function  of  only  the  one  variable,  x.  Now  if  in  addition 
the  moving  load,  P,  is  introduced  with  the  boundary  conditions,  the 
following  differential  equations  are  obtained: 


For  zones  1  and  k 


For  zones  2  and  3 


n  d  w  j.  „u  2  d  w 

D7+PHV     ~2  Sqo 
dx  dx 


(ka) 


k  2 

™  d  w  2  d  w   _,   dw  ,  _        / ,  ,  v 

D     +  pHv  —  -  Cv  —  +  KW  =  qo  (kb) 

dx  dx 


The  boundary  conditions  for  the  indicated  zones  may  be  summarized  as 
follows:   (see  Appendix  I  for  nomenclature) 

MjC-b)  =  0 

Vl(-b)  -  V  (-b) 

w.(-c)  »  0 

w2(-c)  -  0 

i        i 
wx(-c)  -  w2(-c) 

II  II 

w^-c)  -  w2(-c) 

in        in 
w1(-c)  =  w2(-c) 


(a 
(b 
(c 
(d 
(• 
(f 
(8 
(h 
(t 
(J 
(k 


(5) 


w2(0)  -  w3(0) 

w2(0)  =  w^O) 

ii        ii 
w2(0)  =  w3(0) 

w'3"(0)  -  w2(0)  =  P/D 
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(1)  w3(d)  =  o 

(m)  w^(d)  -  0 
(n)  w3(d)  -=  w£(d) 
(o)  w3*(d)  -  .J(d) 

(P)  w3"(d)  -w;'(d) 

(q)  Mu(a-b)  =  0 

(r)  V^(a-b)  =  V  (a-b) 

Solution  of  Problem 
For  the  case  where  A  >0  (Equation  1),  and  the  warped  slab  is  as 
represented  by  Figure  1,  the  deflections  and  stresses  in  the  slab  are 
as  follows:   (see  Appendix  II  for  derivation) 
For  -b  <  x  <   -c 


O  /X,^ 


w  =  A  +  Ax  +  A  cos  rax  +  A,  sin  mx  +  —  (— ) 


'6P 


H 


£  A  (cos  mb  -  cos  mx)  -  A,  (sin  mb  +  sin  mx)~) 


(6) 
(7) 


For  -c  <  x  <  0 
-a.x 


a,x 


°2  = 


6d 

H2 


-axx 


w  =  q  /K  +  e    (B.  sin  b  x  +  B  cos  b  x)  +  e   (B_  sin  b  x  + 
B.  cos  b.x) 

r  i 

B1[  (a2  -  b2,)   sin  b^  -  23^  cos  b2x]  + 

^[^a  sin  b2X  +  (al  "  hP  COS  b2X] 

B-  L(ai    "   bi  )    sin  ^ix  +  ^aibi    cos  blx3  + 
\[(al   -  b2)   cos  blX  -  Pa^   sin  b^]!  +    p(l-Hi)  ^j 


(8) 


+  e 


alX 


(9) 


For  0    <    x   <    d 
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-a,x 


a.x 


w_  =  qQ/K  +  e  (Cx   sin  b^x  +  C     cos  b  x)   +  e        (C_   sin  b  x  + 


C.    cos  b.x) 


(10) 


-6d 

3  2 

°         H 


"alX 


2       .2 


ClC(al   "   V    SlD  V  "   2alb2  C°S  b2XJ  + 


a,x 


C2^2alb2  Sln  V  +  (al   "   b2J   C°S  VJ    +  C   X       C3^al   '   bl*    SiQ  V 


2   ,2 


+  2a  b  cos  b  x]  +  ci,[(ai  "  bi)  cos  bix  "  2a.b  sin  b  xj 


P(l+n) 


H 


(11) 


For  d  <  x  <  (a-b) 

q  2 

Wl+  =  °1   +  °2X  +  °3  C°S  **  +  \   Sin  ""  +  2D   ^ 


(12) 


■69  v 


2  r 


D   £cos  m(a-b)    -   cos  mx  "]+  D,Tsin  m(a-b)    -   sin  bx] 


(13) 


Applying  the  conditions  listed  in  Equation  5>  to  Equations  6,   8, 
10  and  12,  the  following  set  of  non-linear  equations  are  obtained. 


O  /C,^ 


ai  ■  A2° + a3  cos  mc "  kk sin  mc  =  ■  2d  y 


(14) 


a.c 


-a.c. 


e    "   (B     cos  be   -  B.    sin  be)   +  e  (B.    cos  be   -   B_   sin  be) 


"  "   %/K 


(15) 


IT 


a,c 


e       CBi(b2  cos  b2C  +  al   sin  b2C^    "  B2^31   C°S  b2°   "  b2  Sln  b2°^ 

-a   c 
e         [B_(b     cos  be   -  a      sin  be)   +  B.  (a     cos  be  +  b     sin  be)]- 


V 


l  -  m(A  sin  me  +  A,  cos  mc)  ■  -  — r— 

ml) 


(16) 


•  c 


[(a*  -  b2)  cos  b2c  -  2atb2  sin  bgc  ]  -  B]L  [  (aj  -  b|)  sin  bgC  + 


2axb2  cos  b2c] 


+  e 


•v 


B,[2a  bx  cos  be  -   (a2  -  b2)   sin  b  c  ]  + 


o         p  p 

B.    £(a-,    -  b.  )   cos  b.c  +  2a  b.    sin  b.c]    +  m  (A     cos  mc   -  A,    sin  mc) 


o 


Vr 


B2  [i&^l  ~  a  I)   cos  b2c   -  (b|  -   3b2a^)   sin  b2c] 

Bi^(3aib2 '  ai)  sin  V  +  (b2  "  SVl*  cos  V^ 

B3  E(3albl  -  aj)  sin  bxc  -  (bj  -  Sb^)  cos  b^]- 
[(b^  -  3bl3^)  sin  btc  +  (38^  -  aj)  cos  t^c]  + 


-axc 


m  (A     sin  mc  +  A,    cos  mc)   =  0 


B2  +  Brc2-Cr° 


b2(B1  -  Cx)  +  a]L(C2  -  Ck  -   B2  +  Bh)   +  b1(B3  "  C3)  -  0 
2alb2  (C1  "  Bl)  +  (ai  "  b2}  (B2  "  C2)  +   ^l*!  (B3  "  °3)  + 


(17) 


(18) 


(19) 


(20) 


(aj  -  bj)  (B^  -  C^)  -  0 


(21) 
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(b2  "   ^l*    (B1   "  Cl}   "    (3alb2   "  al)    (B2  "  C2}   + 
(bj   -   3bia2)    (b3  -  C3)   +  Oa^  -  .3)    (^   .  C]+) 


P/D 


(22) 


•a.d 


a.d 


e      1   (C,    sin  b  cl  +  C     cos  b  d)  +  e   1   (c„   sin  b,d  +  C,    cos  b,d) 
v   1  2  2  2  3  1  *  1 

-  -  qo/K  (23) 

q         ,2 
D     +  D  d  +  D     cos  md  +  D.    sin  md  =  -°-  (-)  (2*0 


-a,d 


e     l     [C1(b2  cos  b2d  -  ax   sin  bgd)    -  C^  cos  bgd  +  bg  sin  bgd)]     + 


a.d 


e        C  C_(a     sin  b  d  +  b     cos  b  d)   +  C.(a     cos  b  d   -   b.    sin  b  d)] 


qod 
V     +  m(D_   sin  md   -  D,    cos  md)   =  — p- 

vD 


(25) 


m 


-v 


Cl  ^al   "  b2^    Sin  V  '   2alb2  C°S  h2d^+  C2^2alb2  Sin  V  + 


(a!   -  t>2)   cos  b2d] 


+  e 


ad 


C     [(a2  -   b2)   sin  b^  +  2a],b1   cos  b  d] 


+  C^[(a2  -  b2)   cos  bjd  -  2ajb     sin  b  d] 


+  m  (D_  cos  md  + 


D.    sin  md)   =  — — 


o 
mT> 


(26) 


■Vr 


Cl^(3alb2  ■"  al)   8in  V  "    (b2  '   ^2^   °OS  V]  + 


C2^(3alb2  "  al}   ?°S  V  +  (b2  "   3Vl5    SiD  V^ 


+  e 


Vr 


C,[(bJ 


3bia2)   sin  bxd  -   (3axb2  -  a3)   cos  b^]  -  e^£  (^b2  -  a^)   sin  b^  + 


(bl  '  3*1^)  cos  bid] 


-  m  (D-   sin  md  -  D,    cos  md)   ■  0 


(27) 


where  A  =  — 
ra 


V(-b)  sin(mb)  .  /  qo     .  ,,.  .  £T     ,    ,s 

*  in ' +  "X    P(1+^)  nl  cos(mb) 

I  m  D 


mD 


ra 


^'7 

m 


V(-b)  cos  (mb) 

mD 


f<L 


.1X1  I) 


II 


P(l-Hx)  ^  sin(mb) 


^^  sin  m(a-b)  -  |  2fi-  +  p  (l4|i}  £LJ  cos  m(a.b)l 
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(28a) 


(28b) 


(28c) 


and  D,  c  — r 
k         2 
m 


V(a-b) 
mD 


£T 


cos  m(a-b)  +  I  -2-  +  e(l+u)^M  sin  m(a-b)|  (28d) 
1  D 


m 


To  solve  for  the  constants  in  Equations  lA  through  27,  an  initial 
estimate  of  the  solution  was  made  using  the  fully  supported  slab,  then 
a  method  of  functional  interation  equivalent  to  an  N-dimensional  Newton's 
Method  was  used  with  an  IBM  7090  computer  (see  Appendix  III  for  programs). 
With  the  constants  determined,  the  stresses  were  then  obtained  from 
Equations  7,  9,  H  ana"  13«  Using  this  same  procedure,  the  other  deflection 
patterns  resulting  from  the  moving  load  were  analyzed  (see  Figure  2)  and 
the  stresses  were  determined. 


Results 

Due  to  the  large  number  of  variables  involved,  it  is  impractical  to 
present  the  results  for  all  considered  cases.  Instead,  numerical  results 
are  given  for  a  few  typical  cases  to  illustrate  general  trends  and  for 
purposes  of  comparison  with  results  previously  obtained  by  others. 
Typical  programs  are  available  in  Appendix  III  and  they  may  be  used  (or 
modified)  to  obtain  other  values  desired. 

Using  the  left  edge  of  the  slab,  as  origin,  the  moving  load  was 
considered  at  intervals  of  h   ft.  and  stresses  and  deflections  were 
determined  for  conbinations  of  the  following  data. 
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Figure  2.   Typical  Deflection  Patterns 
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u  -  0.15 

6 
E  ■  k   x  10  psi. 

a  =  2^0  ins. 

H  ■  8,  10,  12  Ins. 

K  »  100,  200  pel. 

C  ■  1.5,  2.0 
r    ** 

v  =  0,  20,  U0,  60,  80  mph. 

P  =  125  lbs/in.   (The  value  of  P  was  determined  using  an  axle  load  of 

18,000  lbs.  over  a  pavement  12  ft.  wide) 

q  ■  150.9  pcf. 

tfT-  20,  30,  1*0  deg.  F. 

P  =  6  x  10   inches  per  inch  per  deg.  F. 

Typical  curves  for  deflections  and  stresses  are  presented  in 

Figures  3  through  6  for  three  positions  of  the  load.   For  example,  in 

Figure  3a  the  top  curve  defines  the  deflection  pattern  when  the  moving 

load,  represented  by  an  arrow,  is  at  the  left  edge  of  the  slab.  The 

middle  and  bottom  curves  define  the  deflection  pattern  when  the  load 

is  16  and  32  feet  respectively  from  the  left  edge  of  the  slab.   In 

developing  these  curves,  the  following  data  was  used. 

H  -  8  ins. 

K  -  100,  200  pel. 

C  -  1.5,  2.0 
r       ' 

v  a     kO   mph. 
AT  ■  30  deg.  F. 
Figures  7  and  8  show  how  maximum  deflection  and  maximum  stress  vary 
with  velocity  when  the  difference  between  slab  surfaces  is  20  and  kO 
degrees  Fahrenheit.  A  record  of  the  movement  experienced  by  the  ends 
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of  the  slab  is  presented  in  Figures  IV-1,  IV-2  and  IV-3  (see  Appendix  IV ) 
as  a  function  of  the  load's  position. 

Discussion  of  Part  I* 

Assumptions 

In  developing  the  theory  for  this  part  of  the  analysis,  it  was 
assumed  that  a  pavement  could  be  represented  by  an  array  of  rectangular 
slabs.   Furthermore,  it  was  assumed  that  the  bending  moment  between  slabs 
could  be  neglected  while. a  shearing  force  equivalent  to  that  in  an 
infinite  slab  could  be  used  to  provide  shear  transfer.   This  seemed 
justifiable  because  of  the  relatively  short  depth  of  dowel  embedment 
and  the  general  nature  of  expansion  and  contraction  joints  are  such  as 
to  enable  the  transfer  of  only  very  limited  bending.   For  most  highway 
work,  dowel  bars  are  only  approximately  2  feet  long  compared  to  the 
length  of  a  slab  which  averages  about  kO   feet;  moreover,  they  are, 
generally,  smooth  and  lubricated  at  one  end  in  order  to  maintain  freedom 
of  horizontal  movement  between  slabs.   Under  repetitive  loading,  these 
joints  become  looser  and  conceivably  act  more  like  a  hinge  thus  offering 
little  or  no  moment  transfer.   On  the  other  hand,  substantial  shear 
transfer  could  be  experienced  if  the  joint  opening  is  small  and  the 
deflections  are  large  enough  to  cause  the  development  of  a  bearing 
pressure  between  dowel  and  concrete. 

In  any  case,  it  should  be  noted  that  the  magnitude  of  the  moment 
and  shear  is  really  only  significant  in  the  near  vicinity  of  the  load 


*  The  plan  is  to  treat  each'  part  separately  and  then  summarize  them. 
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itself  (see  Figure  IV-k) .  As  a  result,  the  solution  is  seen  to  be 
influenced  by  the  values  used  for  moment  and  shear  transfer  only  when 
the  load  approaches  the  edges  of  the  slab. 

Results 

Figures  3  and  k   demonstrate  that  for  the  case  studied  here,  the 
damping  ratio,  C  ,  does  not  greatly  influence  the  values  obtained  for 
deflection  and  stress  at  low  to  moderate  velocities;  however,  the  higher 
value  of  C  does  result  in  a  wider  and  less  deflected  trough.*  In  general, 
these  figures  tend  to  indicate  that  the  pattern  of  the  deflection  and 
stress  curves  is  mainly  determined  by  temperature  difference  between 
slab  surfaces  and  by  the  position  of  the  moving  load. 

In  the  case  considered  here,  temperature  differences  cause  the  ends 
of  the  slab  to  curl  upwards  and  become  unsupported,  while  points  midway 
between  the  mid-point  and  the  ends  of  the  slab  experience  an  increase 
in  positive  deflection.   For  the  positions  of  load  shown  in  Figures  3 
through  6,  there  is  an  increase  in  positive  deflection  and  a  decrease  in 
tensile  stress  in  the  near  vicinity  of  the  load;  however,  as  Figures  IV- 1, 
IV-2j and  IV-3,  in  Appendix  IV  indicate,  the  radius  of  influence  (wave 
length  of  deflected  surface)  is  not  only  a  function  of  load  position  but 
also  one  of  velocity.  At  low  velocities,  the  radius  of  influence  is 
small,  but  as  the  velocity  increases,  this  radius  increases  significantly 
behind  the  moving  load.  As  a  matter  of  fact,  this  characteristic  which 
was  also  observed  by  Thompson,  plays  an  important  role  in  the  explanation 
of  Figures  7  a"d  8. 


*  The  depression  under  the  load  is  wider  but  of  smaller  amplitude. 
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For  the  case  of  an  overdamped  pavement  (C  >  1«0),  it  was  initially 
anticipated  that  both  the  maximum  deflection  and  stress  would  decrease 
with  increasing  velocity;  however,  Figures  7  and  8  indicate  that  if 
curling  is  in  evidence  this  may  not  be  the  case.  As  was  stated  earlier, 
the  maximum  positive  deflection  in  an  unloaded  slab  subjected  to  moisture 
and/or  temperature  gradients  which  cause  upward  curling  at  the  ends, 
occurs  somewhere  near  the  l/k   and  3/k   points  of  the  slab.  When  a  moving 
load  is  introduced,  a  wave  train  is  set  up  and  the  deflected  trough 
which  lags  the  load  becomes  wider  as  velocity  increases.   In  other  words, 
the  influence  of  the  deflected  trough  behind  the  moving  load  increases 
with  increasing  velocity. 

Therefore,  if  a  slab  which  is  initially  curled  upwards  at  the  ends 
is  subjected  to  a  moving  load,  there  will  be  a  tendency  for  the  portion 
of  the  slab  behind  the  load  to  become  flatter  and  more  fully  supported 
as  velocity  increases,  see  Figures  IV-1,  IV-2,  IV-3»  As  this  flattening 
occurs,  the  maximum  positive  (downward)  deflection  behind  the  load 
increases  until  the  velocity  reaches  a  value  which  produces  a  reasonably 
flat,  fully  supported  slab;  then,  a  decrease  in  maximum  positive  deflection 
is  experienced  with  further  increases  in  velocity.   Thus  in  a  slab  which 
is  curled  upwards  at  the  ends  the  maximum  positive  deflection  does  not 
occur  at  velocity  =  0,  as  is  the  case  for  fully  supported  slabs  not 
subjected  to  moisture  and/or  temperature  gradients,  but  at  some  velocity 
greater  than  zero. 

As  far  as  stresses  are  concerned,  a  line  of  reasoning  similar  to 
that  used  before  to  explain  the  deflection  pattern,  may  be  employed  to 
account  for  the  increase  in  maximum  stress  with  increasing  velocity  as 
shown  in  Figures  7  and  8.  Increased  velocities  result  in  a  greater 
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tendency  to  flatten  the  slab  and  this  in  turn  causes  an  increase  in  the 
tensile  stress  at  the  top  of  the  slab. 

As  may  be  expected,  the  pattern  as  well  as  the  magnitude  of 
deflections  and  stresses  obtained  depend  on  the  stiffness  of  the  slab, 
the  firmness  of  the  subgrade  material  and  the  temperature  difference 
existing  between  slab  surfaces.   For  the  range  of  velocity  studied, 
maximum  positive  deflection  decreased  as  the  value  of  subgrade  reaction 
increased,  and  increases  were  experienced  as  pavement  thickness  and 
temperature  gradients  became  greater.  However,  as  Figures  7  and  8 
reveal,  the  greater  the  resistance  to  flattening,  the  higher  is  the 
velocity  at  which  the  curled  surface  becomes  flat  enough  to  result  in  a 
decrease  in  maximum  positive  deflection  with  increasing  velocity.   In 
the  case  of  stresses,  the  increase  in  temperature  gradient  from  20  degrees 
to  kO   degrees  Fahrenheit  produced  almost  a  50  percent  increase  in  stress, 
while  the  variation  in  pavement  thickness  and  subgrade  reaction  did  not 
appear  to  have  much  influence.   The  influence  of  the  higher  value  of  C 
on  deflection  and  stress  is  also  small,  however  the  fact  that  it  produces 
a  wider  and  less  deflected  trough  is  quite  evident. 
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PART  II 

In  this  part  of  the  thesis,  the  main  objective  is  to  determine 
what  effect  a  reduction  of  the  subgrade  reaction  would  have  on  the 
stresses  and  deflections  in  a  pavement  subject  to  imposed  vehicular 
loadings.  For  this  case  there  is  little  loss  in  generality  by  assuming 
a  particular  value  of  shear  and  moment  transfer  between  slabs.  The 
procedure  to  be  followed  can  be  adapted  to  any  degree  of  transfer,  hence 
for  the  present  purpose  wherein  there  is  primary  concern  only  for  the 
region  evidencing  subgrade  reduction,  it  is  expedient  to  consider 
primarily  the  case  of  an  infinite  slab.  As  a  special  case,  the  condition 
that  no  shear  and  moment  transfer  exist  between  two  semi-infinite  slabs 
is  also  examined. 

The  pavement  section  for  the  case  of  the  infinite  slab  is  illustrated 
in  Figure  9s-  Here,  the  subgrade  in  its  original  form  is  represented  by 
zones  1  and  3>  while  the  area  over  which  the  reduction  of  subgrade 
reaction  occurs  is  represented  by  zone  2.  The  differential  equation 

describing  the  surface  of  the  pavement  is  given  by  Equation  2,  in  which 

Sw 
p(x. ,  y1 ,  t)  ■  C  sr-  +  Kw.   Using  the  transformation  x  =  (x.  -  vt),  the 

governing  differential  equation  may  be  written  as 

k  2 

_  d  w  .    _  2  d  w   _   dw  ,  „  /„^\ 

D  — jj;  +   pHv  —  -  Cv  —  +  Kw  =  qQ  (29) 

dx  dx 

To  introduce  the  equivalence  of  a  moving  load,  another  zone  must 
be  added  to  the  pavement,  for  example  in  Figure  $b,    zone  1  is  seen  to  be 
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divided  into  tw»  zones,  la  and  lb;  and  the  force  P  is  accounted  for  in 
the  boundary  conditions. 

For  the  case  when  the  moving  load  is  over  zone  1  and  approaching 
zone  2,  Equation  29  is  applied  to  each  of  zones  la,  lb,  2,  and  3,  and 
solutions  are  obtained  using  the  appropriate  subgrade  properties.  The 
boundary  conditions  used  in  evaluating  the  constants  are  as  follows: 
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kl 


Solution  of  Problem 

For  the  case  where  A, ,  A  and  A_  (see  Equation  1)  are  all  greater 
than  zero  and  the  loading  is  as  illustrated  in  Figure  9b,    the  deflections 
and  stresses  in  the  slab  are  as  follows: 

For     x   <    0 


a„x 


w,      =  e       (A_  sin  b,x  +  kc  cos  b,x)  +  qo/K 
la  v    5  16  1   '  o 


(3D 


-6De 
'la"      H2 


31? 


2       .2 


A     ^(a-,    -  b. )   sin  b.x  +  2a..  b.    cos  b  x]  + 


2       ,.2 


A/-  f(a.    -  b. )   cos  b.x  -  2a. b.    sin  b.x  ] 


(32) 


For     0  <    x  <    L. 


a.x 


-a.x 


w       =  e        (B,.   sin  b  x  +  B^  cos  b  x)   +  e  (B„  sin  b  x  +  Bfi  cos  b  x) 


+  qo/K 


(33) 


lb  H2 


r  a.x 


e  l       B     [(a2  -  b2)   sin  b^  +  2a  b     cos  b^  ]   + 


.2       .2 


B6  C  (ai  "  bV  cos  blx  "  2albl  sin  blxJ     + 


"alX 


2       ,_2 


B?[(a^  -  b^)   sin  b2x  -  Sa^  cos  b2x]    + 


2       t2 


B8^al   "  V   COS  b2X  +  2alb2  Sin  b2X^ 


(3»0 


k2 


For  L,  <  x  <  L 

V 


a2x 


w2  "  e  '"  (C5  sin  7jX  +  Cg  cos  7jX)  +  e    (C  sin  y^.  +  CQ  cos  y^s)   + 


i/V^ 


(35) 


-6d 

T     =   — — 

2     2 
H 


a_x 


'  o"    \       i-   p    p 
e  '   K  [  (a2  -  7L)  sin  7^  +  23^  cos  7^]  + 


C6[(a2  -   7j)   cos  7xx  -   28^  sin  7]*]    + 


"V   L    r,J>         2 


c7L(a2  •   ^    Sin  72X  "   2a272  C°S  72X^  + 


Cg  [2a272  sin  72x  +  (ag  -  7g)   cos  7gx] 


(36) 


For  x    >  L2 


■a,x 


w     =  e  (D„  sin  b  x  +  Dg  cos  b  x)   +  qo/K 


(37) 


ztt,     -a,x 

■  -oD  .      1 

*       H 


2       ,2 


D?[  (a J  -  b2)    sin  b2x  -   2aib2  cos  bgx  ]    + 


2       ,2 


(Sl   "  b2^    °OS  b2X  +  2*)}>2  Sln  b2XJ 


(38) 


Applying  the  conditions  listed  in  Equation  30  to  Equations  31,  33* 
35  and  37,  the  following  set  of  linear  equations  are  obtained. 


A6  "  B6  "  B8  =  ° 


(39) 


^3 


blA5  +  alA6  "  blB5  '  alB6  "  b2B7  +  alB8  "  °  (^) 

2alblA5  +  <al  "  bl>  A6  "  2alblB5  "  (al  "  bl>  B6  * 

(a2  -  b2)  Bg  +  23^^  =  0  (1H) 

(bj  -  3bla2)  A5  -  (a3  -  ^b2)  A6  -  (bj  -  35^)  ^   + 

(al  "  3albl)  B6  +  (^2  "  b2)  B7  +  (3alb2  "  al}  B8  =  P/D        (te) 

a  L  -a  L 

e  (Br   sin  b.L.    +  B^   cos  b.L.  )   +  e  a(b_  sin  bl.    +  BQ  cos  b  L   )    - 

x5  116  1   1'  v7  218  2  1' 

32L1  'a2Ll 

e  '      (C5  sin  7^  +  C^  cos  7^)   -  e      '     (C?  sin  7^  +  CQ  cos  7^) 


alLl 

e 


[|  B[-(b1    cos  b.L     +  a      sin  b.L.)   +  B/-(a..    cos  b.L. 


"alLl 

bx   sin  ^1^)]  +  e  L^b2  cos  b2Ll   "  ai   sin  Vl^ 


BQ(b2  sin  b2Lx  +  3j   cos  b^)]"   e   ''   \cAv1   cos   7^ 


+  &2  sin  7XLX)  +  C6(a2  cos  7^   -  y1   sin  7^)  J+  e  [C8^2  Sin  72L1 


+  a2  cos  7^)   -  C7(72  cos  7^   -  a2  sin  72\)1:=  0  (hh) 


alLl 


B5  C(ai  "  bi)   8in  b1L1  +  2a-^x  cos     biLi  ]  + 


2   t2 


Bg  [(a^  -  bp  cos  b^  -  2a1b1  sin  b^] 


"alLl 


B?[(a^  -  b2)  sin  b2Lx  -  23^  cos  1^]  + 


Bg  C(a^  -  b|)  cos  bgLj  +  2axb2  sin  b^  ] 


-   e 


a2Ll 


S  ^a2   "   V    Sin  71L1  +  2Vl   C°S  71L1^+ 


2  2  -J 

C6[(a?  -   7X)   cos  7XLL   -  23^   sin  7^A' 


•a2Ll 


2    2 
C7  t^a2  "  72*    Sin  72L1  "  2a272  C°S  72LJ  + 


1* 


C8^2a272  Sin  72L1  +  ^2  "   72^   C°S  72Ll3 


0 


(*5) 


*1L1 


B5C(a^  -  3^^)  sin  b^  -  (bj  -  frff)   cos  b^  ]  + 


B6[  (b3  -  3bia2)  sin  b^  +  (a3  -  3aib2)  cos  b^J 


"alLl 


B  [(3a^b2  -  b3)  cos  b2Lx  +  (33^  -  a3)  sin  b^  ]  + 


Bg  [(3*^2  -  a3)  cos  b2Lx  +  (b3  -  33^)  sin  b^  ] 


.a2Ll 


2\   ,_   ,    ,3   -,.  -2, 


c5  E(a2  _  ^V  Sin  71L1  '  (71  "  37ia2^  C°S  71L1  J 


+ 


C6  E(71  "  37132)  Sin  71L1  +  (a2  "  3a27l)  COS  71L1  J 


-a  L 
2  1 


c7  C(3a272  "  72)   cos  72LX  +  (3a272  -  a3)   sin  y^J  + 


^5 


C8  C(3a272  "  a3)   cos   7^  +  {y\  -   3a272)    sin  7^] 


0 


(^6) 


a  L  -a  L0 

e  2  2(C5  sin  7XL2  +  Cg  cos  7^)  +  e  2  2(C?  sin  7^  +  Cg  cos  7^) 


e  &1  2(D?  sin  b2L2  +  Dg  cos  b^)  =  %H^) 


(*7) 


h6 


S2L2 


e    [  C5  (jj  cos  7^  +  a2  sin  7^)  +  Cg  (a2  cos  7^  - 


2  2 


x   sin  7^)]  +  e     [  C7  (72  cos  7^  -  a2  sin  7^)  - 


-a,L, 


Cg  (72  sin  72L2  +  ag  cos  7^>2)~\     -e     l  2  L^i^   cos  b^  - 


aY   sin  b2L2)  -  Dg  (bg  sin  b^  +  a^   cos  b^)  ]  =  0 


(1*8) 


a2L2 


C5  ^a2  "  71^  Sln  71L2  +  2a27l  C°S  71L2  ^  + 


C6  [(a|  -  7\)    cos  7lL2  -  2a27l  sin  y^J^ 


+  e 


•a2L2 


.2    2 


C?[(a2  -  7g)  sin  7^ 


-  2a272  cos  72L2]  +  CQ  [2a272  sin  7^  +  (a£  -  7£)  cos  72L2] 


"alL2 


D7  t(al  "  b2J  Sin  b2L2  '  2aib2  C°S  b2L2^  + 


D8C2aib2  Sln  b2L2  +  (ai  "  b2)  C°S  b2L2^ 


0 


(^9) 


vr 


.  2  2 


C5  [(a2  "  3a27l^  Sin  71L2  '  ^71  "  37la2^  C°S  71L2^J  + 


C6^^7l   '   3V2}  sin  ^1L2  +  (a2  '  3a27l)  COS  71L23 


"a2L2 


c7  [(3*2^2   '   72>   C°S   72L2  +  (3a272  "  a2>    sin  72L2]  + 


C8[  (3a272  "  a|)   cos  y^  +  {y\  -   yfy2)   sin  y£>A   - 

DT  ^3a^b2  -  b3)   cos  b2L2  +  (38^  -  a3)   sin  b2L2]  + 


"alL2 


D8C  (3alb2  "  V   COS  b2L2  +  (b2  "   3a?2)    Sin  b2L2^ 


0 


(50) 


The  constants  in  Equations  39  through  50  were  evaluated  using  the 
Method  of  Crout  Reduction  (25)  with  an  IBM  7090  computer.   In  addition, 
a  column  pivot  search  was  performed  during  reduction  in  order  to  improve 
accuracy.  After  evaluating  the  constants,  the  deflections  and  stresses 
were  determined.  The  complete  solution  to  the  problem  was  obtained  by 
applying  this  same  procedure  with  suitable  modifications  to  the  schemes 
shown  in  Figure  9C  and  9d.  The  general  solution  for  the  case  where  A. 
and  A3  >  0  but  ^  <  0  is  given  in  Appendix  II. 

To  examine  the  effects  of  a  loss  in  load  transfer  at  a  joint  or 
crack  assumed  to  exist  at  the  mid-point  of  the  weakened  zone,  the  special 
case , illustrated  in  Figure  10  was  investigated.  Here,  the  moment  and 
shear  acting  at  points  e  and  f,  are  assumed  to  be  respectively  0$  and 
505&  of  that  which  would  exist  in  a  corresponding  section  of  an  infinite 
slab.  Using  the  same  procedure  detailed  above,  a  set  of  linear  equations 
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were  developed.  The  constants  were  evaluated  and  the  stresses  and 
deflections  determined. 

Results 
To  investigate  the  effect  of  a  reduction  in  subgrade  reaction, 
numerical  results  were  obtained  for  combinations  of  the  following  data: 
u  =  0.15 
E  =  k   x  10  psi 
£   =  2k,   48,  96  ins. 
H  =  8,  10,  12  ins. 
K  -  100,  200  pci 


K.=  0,  25,  50,  75,  100,  150,  200  pci 


C  «=  1.5,  2.0 
r     ' 

v  =  0,  20,  40,  60,   80  mph 
P  =  125  lbs/in. 

qQ=  150.9  pcf 

Some  typical  curves  (in  non-dimensional  form)  for  displacements  and 
stresses  for  three  positions  of  the  load  are  presented  in  Figures  11,  12 
and  13. 

In  Figures  11  and  13,  the  moving  load  (represented  by  the  vector, 
P)  is  2  feet  away  from  the  nearest  point  of  the  zone  of  reduced  subgrade 
reaction,  K- ,  while  in  Figure  12,  the  load  is  at  the  mid-point  of  this 
zone. 

Values  for  the  deflection  and  stress  under  a  static  load,  P,  are 
given  in  Table  1.  These  values  may  be  used  in  Figures  11  to  13  with  the 
appropriate  amplification  ratios,  to  obtain  the  magnitudes  of  the  deflection 
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Table  1 
Deflection  and  Stress  Under  a  Static  Load  of  125  lbs/in 
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K     =  100 
o 

pci 

K     =  200 
0 

^ci 

H  ins 

w     ins 
o 

.0191M+ 

a     psi 
0  r 

150.607 

w     ins 
0 

.010722 

a     ins 
0 

126. 6^+5 

8 

10 

.019017 

113.9^8 

.oioWl 

95.819 

12 

.019^9 

90.726 

.010573 

76.291 

5^ 


and  stress  at  a  point.  For  example:  For  a  pavement  thickness  of  8 
inches  and  K  of  200  pci,  the  maximum  stress  under  a  static  load  of  125 
lbs/in,  is  seen  from  Table  1  to  be  126.645  psi.  From  Figure  12,  for  a 
velocity  of  40  mph,  VL  /K  ■  0.5  and  C  =  1.5,  the  stress  amplification 
for  a  point  24  inches  behind  a  load  situated  at  the  mid-point  of  the 
weakened  region,  is  0.431.   The  magnitude  of  the  actual  stress  is  then 
0.431  x  126.645  -  54.64  psi. 

Figures  14  and  15  illustrate  the  variation  of  maximum  stress  and 
deflection  with  velocity,  while  Figures  16,  17,  18  and  19  reveal  the 
influence  of  the  parametric  length,  £  ,     Similar  plots  are  presented  in 
Figures  20  through  23  for  the  parameter,  H. 

Curves  similar  to  those  given  in  Figures  11,  12  and  13  are  presented 
in  Figures  24,  25  and  26  for  the  special  case  where  there  is  only  50$ 
load  transfer  and  no  moment  transfer  across  a  discontinuity  in  the 
pavement's  surface. 
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Figure  14a  Maximum  Deflection  and  Stress 

Amplification  vrs  Velocity  (KQ=  200  pci., 


Gr=1.5) 
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Figure  22a.     Maximum  Deflection  Amplification  versus   Slab 
Thickness   (K  «100  pci.,   C  =1.5) 
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Figure  22b.  Maximum  Stress  Amplification  versus  Slab 
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78 
Discussion  of  Part  II 

As  shown  in  Figures  11,  12  and  13,  there  is  an  increase  in 

positive  deflection  amplification  as  the  value  of  K, /K  decreases.  For 

1  o 

the  load  positions  shown,  the  increase  is  greatest  when  the  load  is  over 

the  zone  of  reduced  subgrade  reaction  and  least  when  the  load  is  passed 

this  zone.  For  example:  for  a  value  of  K./K  =  0.5,  the  load  at  raid-point 

(Figure  12)  yields  a  deflection  amplification  ratio  ss  1.28,  the  load  2  ft, 

behind  the  zone  of  reduced-K  (Figure  11)  gives  approximately  1.1  and  the 

load  2  ft.  in  front  of  this  zone  yields  ^  1.07  (Figure  13). 

As  may  be  expected,  when  the  load  is  moving  over  a  homogeneous 

subgrade  material,  the  point  of  maximum  positive  deflection  lags  the 

position  of  the  load;  however,  if  there  is  a  zone  of  reduced  subgrade 

reaction  and  the  load  is  approaching  this  zone,  as  is  shown  in  Figure  11, 

the  maximum  deflection  can  lead  the  position  of  the  load  at  low  values 

of  K,  /K  .  For  the  cases  when  the  load  is  over  the  weakened  region  or  has 
1  o  ° 

passed  it,  maximum  positive  deflection  lags  the  load.  This  characteristic 
tends  to  increase  as  the  load  moves  out  of  the  weakened  area. 

As  regards  to  stresses,  a  reduction  in  the  value  of  K  does  not 
appear  to  influence  the  maximum  stress  very  much  unless  the  load  is 
within  the  area  of  reduced-K;  however,  as  Figures  11  and  13  demonstrate, 
it  can  effect  significantly  the  shape  of  the  curve  at  very  low  values  of 

w 

For  the  3  positions  of  load  considered,  the  maximum  stress,  which 
always  occurred  under  the  moving  load,  reached  its  greatest  value  when 
the  load  was  over  the  weakened  zone;  here,  stress  increased  as  K. /K 
decreased.  In  the  case  when  the  load  is  approaching  the  weakened  zone 
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stress  is  again  seen  to  increase  with  decreasing  ratio,  KL./K  ,  but  after 
the  load  has  passed  this  zone,  the  trend  is  reversed  and  the  lower  the 
value  of  K-/K  ,  the  lower  the  stress  under  the  load.  This  means  that  if 
a  slab  is  lying  on  a  subgrade  which  exhibits  a  soft  spot,  the  stress 
experienced  by  the  slab  is  greater  when  the  load  is  moving  towards  the 
center  of  this  weakened  zone  than  when  it  (the  load)  is  moving  away. 
Consequently,  if  failure  due  to  overstressing  does  occur,  the  signs  of 
distress  should  appear  somewhere  between  the  center  and  the  "back"  edge 
of  the  soft  zone. 

Figures  lU  and  15  demonstrate  that  except  for  very  low  values  of 
K-/K  ,  deflection  and  stress  amplification  decreases  with  increasing 
velocity,  however  the  decrease  appears  to  be  more  pronounced  for  the 
smaller  value  of  initial  subgrade  reaction,  K  .  Here,  as  in  Part  II,  the 
higher  value  of  C  yields  the  lower  value  of  maximum  deflection  and  stress. 

The  influence  of  the  length  of  weakened  region  is  illustrated  in 
Figures  16,  17,  18,  and  19.   In  these  figures,  both  the  maximum  stress 
and  deflection  are  seen  to  increase  for  all  values  of  K./K  less  than  1, 
as  the  region  of  reduced  subgrade  reaction  becomes  larger.   The  stiff er 
the  subgrade  material  is  initially,  the  greater  the  increase  in  deflection 
and  stress.   However,  in  the  case  of  stress,  this  is  only  discernible 
at  low  values  of  K-./K  . 

In  general,  the  influence  of  the  thickness  of  the  pavement  upon 
stress  and  deflection  amplification  was  only  significant  at  low  ratios 
of  1C./K  .  As  is  shown  in  Figures  20  through  23,  deflection  amplification 
decreases  as  thickness  increases,  when  the  value  of  K-/K  is  approximately 
equal  to  0.7  or  less;  while  stress  amplification  remains  relatively 
insensitive  to  variation  in  thickness.  Again,  these  plots  show  that  for 


3o 

constant  K  /K  ratios  less  than  1,  the  initially  stiffer  gubgrade  leads 
to  the  greater  increase  in  stress  and  deflection.  This  suggests  that  if 
there  is  the  possibility  for  the  development  of  soft  spots  in  the  rubgrade 
material,  it  may  be  better  to  avoid  the  use  of  stiff  subgrades. 

In  Figures  2h,    25,  and  26,  there  is  a  clear  indication  of  what 
happen,  when,  in  addition  to  the  reduction  in  subgrade  reaction,  50^  of 
load  transfer  is  lost.  Deflections  are  substantially  increased  near  the 
point  of  discontinuity,  expecially  when  the  load  is  over  the  weakened 
region.  For  instance,  in  Figure  25a,  for  YL  /K  -  0.5,  the  deflection 
anpjif ication  for  a  point  under  the  ooving  load  is  >.8>  compared  to  1.26 
in  Figure  12a.  Another  important  aspect  to  note,  i3  the  large  relative 
movement  which  occurs  between  the  edges  of  the  slab  (i.e.  edges  situated  at 
the  mid-point  of  the  weakened  region).  This  -.ovement  not  only  causes 
bumpy  driving,  but  nay  also  lead  to  further  pavement  distress. 

As  for  stresses,  Figure  15b  demonstrates  that  for  the  special 
case  considered,  the  maximum  stress  need  not  occur  under  the  moving  lc 
In  fact,  the  maximum  stress  experienced  in  this  case  is  not  only  behind 
the  load,  but  is  also  correspondingly  higher  than  that  obtained  in  Figure 
12  -   For  the  other  two  positions  of  the  load  (Figures  2^b  and  26b),  t 
use  of  no-moment -transfer  and  yCffi   shear  transfer  did  not  greatly  influence 
the  values  previously  obtained  in  Figures  11  and  13  tot  wmm   stress,  as 

long  as  the  ratio  of  K_/K  was  high  (approximately  =  0  .     owever,  when 
the  ratio  of  K  /K  was  lowered  to  0.25,  significant  reduction  in  stress 
was  experienced.  Obviously,  at  high  values  of  L/K  ,  the  load  in  t: 
case  has  to  be  fairly  close  to  the  weakened  region  before  conditions 
existing  at  the  mid-point  of  the  region  become  important. 
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SUMMARY 

In  Part  I  of  this  thesis,  it  is  shown  that  when  a  pavement  is 
subjected  to  upward  curling,  it  is  possible  to  experience  an  increase 
in  maximum  positive  deflection  as  velocity  is  increased  and  then  a 
decrease  (in  maximum  positive  deflection)  with  further  increases  of 
velocity.  The  velocity  at  which  the  decrease  in  maximum  positive 
deflection  sets  in  seems  to  depend  on  several  factors  which  include 
the  thickness  of  the  slab,  the  temperature  difference  between  its 
surfaces,  the  stiffness  of  the  foundation  and  the  degree  of  damping. 
In  general,  it  appears  that  higher  velocities  are  needed  to  effect  a 
decrease  in  maximum  positive  deflection  as  slab  thickness,  stiffness  of 
foundation  and  temperature  increase,  and  as  the  degree  of  damping  decreases. 

In  the  case  of  stresses,  maximum  stress  increased  with  increasing 
velocity.  The  thinner  the  pavement  and  the  lower  the  value  of  subgrade 
reaction  and  damping  coefficient,  the  higher  the  stress;  however,  there 
was  not  a  great  difference  in  the  values  of  stress  obtained.  What  seemed 
to  matter  most,  was  the  difference  in  temperature  between  the  surfaces 
of  slab.   Increases  in  the  temperature  difference  resulted  in  large 
increases  in  stress  and  maximum  positive  deflection,  for  all  values  of 
velocity  studied.  This  observation  tends  to  indicate  temperature  difference 
between  slab  surfaces  as  the  over-riding  factor  governing  the  magnitude 
of  stress  and  deflection  which  may  be  obtained,  while  factors  such  as 
velocity,  thickness  of  pavement,  modulus  of  subgrade  reaction  and  degree 
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of  damping  act  more  or  less  to  exaggerate  or  minimize  the  pavement 
distress  caused  by  temperature  and/or  moisture  gradients. 

In  Part  II,  where  the  influence  of  a  reduction  in  subgrade  reaction 
was  studied,  it  was  clearly  shown  that  both  maximum  deflection  and  maximum 
stress  increased  as  the  region  became  weaker.  The  increase  experienced 
was  quite  pronounced  when  the  load  was  within  the  weakened  region  and 
moving  towards  the  center.   Thus,  if  a  pavement  is  designed  with  a 
particular  value  of  subgrade  reaction,  and  a  weakened  zone  develops  due 
to  the  softening  of  ground  (as  may  be  the  case  during  Spring  thaw)  the 
stresses  produced  in  the  slab  in  the  near  vicinity  of  the  weakened  zone 
will  very  likely  be  higher  than  anticipated. 

If  there  is  significant  reduction  of  subgrade  support  (even  to 
within  75$  or"  tne  original  value  of  K  ),  maximum  deflection  and  stress 
should  still  follow  the  well  known  trend  and  decrease  with  increasing 
velocity;  however,  when  there  is  complete  loss  of  subgrade  support,  both 
the  maximum  stress  and  deflection  can  increase  with  increasing  velocity. 
For  the  hypothetical  case  where  pavements  of  equal  thickness  are  built 
on  different  subgrade  material,  it  appears  that  for  a  constant  percentage 
loss  of  subgrade  support,  the  pavement  built  on  the  stiffer  subgrade 
material  should  experience  a  greater  increase  in  stress  and  deflection. 
In  the  case  where  there  is  also  a  loss  of  load  transfer,  deflections  as 
well  as  stresses  may  be  substantially  increased  and  large  relative  move- 
ment may  be  experienced  near  the  points  of  discontinuity  in  the  surface 
of  the  slab. 
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CONCLUSIONS 

1.  On  the  basis  of  the  assumptions  stated  herein,  a  theory  has 
been  developed  whereby  stresses  and  deflections  can  be  computed  for 
finite  rectangular  slabs  subject  to  ambient  and  imposed  moving  loads. 
The  base  may  exhibit  time  dependent  effects  and  may  also  be  non- 
homogeneous  . 

2.  Contrary  to  common  opinion,  the  theory  demonstrated  that  an 
increase  in  velocity  can  produce  an  increase  in  deflection  and  stress. 

3-  Of  all  the  variables  considered  herein  (temperature  differences, 
velocities,  thicknesses  of  pavement,  moduli  of  subgrade  reaction  and 
degrees  of  damping),  the  temperature  difference  between  slab  surfaces 
was  shown  to  be  the  over-riding  factor  governing  the  magnitude  of 
maximum  stress  and  deflection. 

h.     There  is  reason  to  believe  that  any  meaningful  interpretation 
of  the  performance  of  pavements  as  determined  by  measured  strains  and/or 
deflections  of  the  slab,  must  give  due  regard  to  the  effects  of  warping. 

5-  The  results  of  this  study  indicate  that  stresses  and  deflections 
within  a  slab  can  be  sensitive  to  localized  reductions  in  subgrade  support. 
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SUGGESTIONS  FOR  FURTHER  RESEARCH 

The  procedure  developed  here  could  equally  be  applied  to  the  case 
where  the  temperature  at  the  top  of  the  slab  is  greater  than  that  at 
the  bottom.   Furthermore,  other  models  such  as  the  Maxwell  and  the 
Standard  Solid  may  also  be  used  to  simulate  support  conditions. 

At  the  present  time,  there  is  very  little  known  on  the  range  of 
values  to  be  used  in  viscoelastic  models  and  therefore  a  full  scaled 
experiment  is  in  order.   Such  an  experiment  should  not  only  measure 
stresses  and  deflections,  but  should  in  particular  pay  attention  to  the 
prevailing  boundary  conditions.   For  example,  some  attempt  should  be 
made  to  evaluate  the  conditions  at  the  joints  in  highway  slabs. 

The  problem  of  a  slab  lying  on  a  viscoelastic  material  and  subject 
to  warping  as  well  as  to  a  series  of  moving  concentrated  loads  is  also 
worthy  of  future  study. 
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Nomenclature 

AT  ■  difference  in  temperature  between  the  surfaces  of  the  slab 

t  =  time 

w  =  mid-plane  deflection  of  the  slab  (positive  downward) 

x- ,y.  ■  fixed  coordinates 

x,  y  a   coordinates  with  respect  to  moving  load 

E  =  Young's  modulus 

u  =  Poisson's  ratio 

H  =  slab  thickness 

q  »  uniform  distributed  load  due  to  the  weight  of  the  slab 

P  -   concentrated  line  load 

p  =  density  of  the  slab 

K  =  modulus  of  subgrade  reaction 

C  =  damping  coefficient 

v      =  velocity  of  moving  load 

EH 
D      =»  —  ■  flexural  rigidity  of  the  slab 

12(l-u2) 

a      ■  length  of  slab 

/      ■  length  of  region  of  reduced  subgrade  reaction 

L.  ,L   «■  distances  from  load  to  edges  of  region  of  reduced  subgrade 
reaction 

b,c,d  ■  distances  from  load  to  edges  of  zones  in  partially  supported 
slab 

* 

C      ■  ratio  of  the  damping  coefficient  to  that  of  the  critical 
damping  coefficient 


V(x)   «  shear  at  any  point  x 
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M(x)  ■  bending  moment  at  any  point  x 

V(x)  ■  shear  at  any  point  x  in  infinite  slab 

a(x)  ■  tensile  stress  at  any  point  x 

p  =  linear  coefficient  of  thermal  expansion 
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APPENDIX  II 

GENERAL  SOLUTION 
Part  I 
The  general  differential  equation  for  the  deflection  of  a  thin 
rectangular  plate  resting  on  a  viscoelastic  foundation  and  subjected  to 
a  moving  line  load,  P,  is  (42): 

D^l  +  PH  ^  +  C  ^  +  Kw  =  qo  +  P  (.  ,t)  (1) 

dx.         St 

Employing  the  transformation  x  =  (x1  -  vt),  and  introducing  the  moving 

load,  P,  with  the  boundary  conditions,  Equation  (l)  becomes: 

k  2 

_  d  w  -,  2  d  w   _  dw  ,  _                       ,_» 

D  — r—  +  pHv   — -  -  Cv  —  +  Kw  =  q                     (2) 

,   k  .  2      dx        Mo                   v  ' 

dx  dx 

The  solution  of  Equation  (2)  as  given  by  Thompson  (59)  is: 

For  A  >  0 

q      "aix  aix 

w  =  -— ■  +  e     (N..sin  b  x  +  N0cos  b  x)  +  e  "   (N  sin  b.x  +  Nj  cos  Kx)  (3a) 

For  A  <  0 

q       a,x      a_,x    "aix 
w  =  —   +  N  e    +  N^e  5     +   e     (N  sin  b  x  +  NQcos  b  x)  (3b) 

2        k 
,      6    .   m       k    ,    ,m  K  n   2   /Cvx2   _  /,,  \ 

where  «x  +  _  &i  +   (^  -  ,_)  ax  -  (^)  =  0  (lw) 


'k-'i+l^l'^-T  ^ 


/  Cv      2   m2 
a3  =  3l  Ih^-   31  "2 


(4c) 


h  -  k  -  & + 4  <*> 
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(fcf) 

For  the  region  where  the  plate  is  unsupported  Equation  (2)  reduces 

to: 

k  2 

r,  d  w  i  »  „  2  d  w  /  c  \ 

D  — £  +PHv  — -  =   qQ  (5) 

dx         dx 


b2-l 

/   2    ,      Cv          m2 
K   +  1^   +  2 

.      2      P  Hv2 
and  m     -  —  ■ 

qo  ,x^2 


The  solution  to  Equation  (5)  is: 

w  -  N  +  N1Qx  +  Nn-  cos  (mx)  +  N^  sin  (rax)  +  -^  (^  (6) 

where  m  is  given  by  Equation  (Kt ) . 

To  obtain  the  solution  for  the  case  described  by  Figure  1,  the 
following  conditions  must  be  satisfied: 

(a)  M^-b)  =  0 

(b)  vx(-b)  =  V(-b) 

(c)  Wl(-c)  =  0 

(d)  w2(-c)  =  0 

(e)  wj(-c)  =  Wg(-c) 

(f)  wj(-c)  -Wg(-C)  (7) 

(g)  w"(-c)  =  w™(-c) 
(h)  w2(0)  -  «3(0) 

(i)  w2(0)  =  w^O) 

(j)  «2(0)  =  w3(0) 

(k)   w'''(0)  -  w'2'(0)  =  P/D 
(1)   w3(d)  =  0 
(m)   w^(d)  =-0 


<* 


(n)  w'3(d)  -w'u(d) 

(o)  w"3(d)  =  v"^(d) 

(p)  w'"3(d)  =  w'«  (d) 

(q)  M^  (a-b)  =  0 

(r)  V^  (a-b)  =  V  (a-b) 

For  the  case  where  A  >  0,  the  deflection  of  the  plate  may  be 

expressed  as  follows: 

qo  x  2 

Wl  =  Al  +  A2X  +  A3COS  (""J  +  Al,.sin  (n")  +  2D  @  ^ 

q     -a.x  a.x 

w  ■  —  +  e     (B.sin  b  x  +  B  cos  b  x)  +  e    (B  sin  b  x  + 

B^cos  b^x)  (9) 


<U    "aTX 


a.x 


w  =  —  +   e     (C  sin  b  x  +  C  cos  b  x)  +  e    (C  sin  b  x  + 

C,  cos  b  x)  (10) 


q    2 

O  /Xn 


w^  =  Dl   +   D2x  +  D3cos  (mx)  +  D^sin  (nx)  +  —  (-)  (ll) 


The  bending  moment  in  a  thin  plate  subject  to  a  gradient  AT/H  (6l)  is: 

& 


p 

m-  -d  m+  p(i+^)^i  <i2> 

Ldx 


For  M1(-b)  =  0,  Equation  (7a) 

p  p  a  AT 

Am  cos  (mb)  -  A.m  sin  (mb)  =  —j-  +  P(l+u)  —  (13) 

ml) 

The  shear  in  a  thin  plate  subject  to  a  gradient  AT/H  (6l)  is: 

v  = -°  3H7I +  ?  (1^>  *t]  <l*> 

^dx 
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As  the  gradient  is  assumed  to  be  independent  of  x,  Equation  (lU)  reduces 
to: 


V  -  -  D  A) 

dxJ 


For  V.(-b)  =  V  (-b),  Equation  (7b) 

A,m  sin  (mb)  +  A, m  cos  (mb)  =  — *rr 


Solving  Equations  (13)  and  (l6) 
1 


A3=  2l 

m  L 


V(-b),ln(mb)  J!»    p(14v)ffi) 

ImT) 


mD 


m 


V(-b)cos(rab) 
mD 


M 


H 


AT 


cos(mb) 


-f-  +  e(l-H0~l  sin(mb)! 


(15) 


(16) 


(17a) 


(17b) 


For  M.  (a-b)  =  0,  Equation  (7q) 


?  2  o  AT 

Dm  cos  m(a-b)  +  D.m  sin  m(a-b)  =  ~x~  +    g(l+u)  TjT- 

^  ml) 


For  V,  (a-b)  -  V(a-b),  Equation  (7r) 


•D_m  sin  m(a-b)  +  D.nr  cos  m(a-b)  =  — *"r 


V(a-b) 


(18) 


(19) 


Solving  Equations  (18)  and  (19) 


m 


V(a-b) 


AT 


mD 


sin  m(a-b)  -  -5-  +  g(l+n)  —  cos  m(a-b) 


ml) 


H 


m 


V(a-b) 


mD 


cos  m(a-b)  + 


'^-+  P(l+u)^|  sinm(a-b) 


(20a) 


(20b) 
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Applying  Equations  (7c)  through  (7p)  to  Equations  (8),  (9),  (10) 
and  (11),  the  following  set  of  non-linear  equations  are  obtained. 


qo     c   2 
Al   "  A2C  +  A3  cos(mc)    "  \  sin(mc)   -   -  —  (-) 


(21) 


a.c 


a,c 


e        (Bp  cos  be   -   B1    sin  be)  +  e  (B.    cos   b.c   -  B_   sin  b.c)   « 


VK 


(22) 


a,c 


1     [B1(b2  cos  b2c  +  a1   sin  b^)    -  B2(a1   cos  b^   -   bg  sin  bgc)] 


-a,c 


+  e  £B_(b     cos   b..c   -  a     sin  b.c)   +  B.  (b     sin  be  +  a      cos   b_c)~| 


A.     -  Am  sin  mc   -  A,  m  cos  mc 


V 


m 


*D 


(23) 


V 


2        ,2 


2        ,2, 


B2[  (aj   -   b2)   cos   b2c   -   23^  sin  b^  ]  -  Bx  [(a^  -   b^)    sin  bgC 


+  2atb2  cos  b2c] 


+  e 


-V 


B,    [(a      -   b.)   cos  be  +  2a   b     sin  be] 


2        ^2 


+  B,[2a  b     cos   be   -   (a      -   b.  )   sin  be]       +  m  (A     cos  mc   -  A,    sin  mc) 


ra 


% 


(2h) 
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V 


B  [Oa^  -  aj)  cos  b2c  -  (b^  -  3a^b2)  sin  bgc  ]■ 


Bx  [(3"j[b2  '  b2)  C0S  V  "  (3alb2  "  31}  Sin  V^ 


-a.c 


1   [b3  [(a3  -  38^)  sin  b^  +  (bj  -  3*^)  cos  b^]* 

\  [(bl  '  ^lV  Sln  bl°  "  (S1  "  ^^  C°S  ^^  +  " 
m3  (A  sin  mc  +  A^  cos  mc)  =  0 


(25) 


B2  +  \   -  C2   '   Ck  =  ° 


b2(B1    -   Cx)   +  al(C2   -   Ch    -   B2  +  B^)    +  h^   -  Cg)   -  Q 
2aib2(Cl   "  Bl)   +  (^   "   b2}    (B2   "  C2)   +  23lbl   U3   ■  °3)   + 


(a*   -  bj)    (BU   -  C^)   -  0 


2       „3> 


(b|  -   3b2a^)    (Bl   -  Cx)    -    (3-jbl   -  .J)    (B2  -  C2)   + 

(b'3  -  3^)  (b3  -  c3)  +  (Sa^  -  aj)   (b^  -  ck)  -  I 

-a  d  aid 

e  *   (Cx  sin  b2d  +  C2  cos  bgd)  +  e    (C3  sin  b^  + 


(26) 


(27) 


(28) 


(29) 


C^  cos  bjd)  -  -  qQ/K 


(30) 
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D  +  D  d  +  D   cos  md  +  D,  sin  md  »  -—  (— ) 


(3D 


-a.d 


e  l     Cc1(b2  cos  V  ■  ai  sln  b2d)  "  C2^al  C°S  V  +  b2  Sin  V^  + 


a.d 


1     CCo^al    sln  bid  +  bi    cos  bid)   +  Ck(ai    cos   bid   "   bi    sin  bid)]  ' 


%d 


D  +  ra  (D   sin  md  -  D,  cos  md)  ■  -^— 

ml) 


(32) 


-a.d  r 


2   ,2 


Cl  C(al  "  V   sin  V  '  2alb2  C°S  V  J+  C2  ^2alb2  Sin  V 


+  (al   "  b2^   cos  b2d] 


+  e 


V 


C    [(a2  -  b2)   sin  b^H-  2a1b1   cos  b  d] 


(a      -  b   )   cos  b  d  -  2a  b     sin  b  d]     +  m  (D     cos  md  +  D,    sin  md) 


no 

ml) 


"aid 


Cl  t(3aib2  "  *V    Sln  V 


(b2   "   ^l*    C°S  b2d^  + 


C2  ^(3aib2  "  31)   COS  V  +  (b2  "   3Vl>    Sin  V^     + 

Cu[(b3  -  3^)   sin  bxd  -   (3aLb2  -  a3)   cos  b^] 


aid 


C3  C(3axb2  -  a3)   sin  b^  +  (b3  -  S^a2)   cos  bjd] 


(33) 


•3 

m  (D-   sin  md   -  D»    cos  md)   =  0 


W 
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With  the  use  of  N- dimensional  Newton's  Method  (27),  the  constants 
in  Equations  21  through  3^  may  be  determined  and  the  deflections  obtained 
from  Equations  8,  9,  10  and  11. 

Using  Equation  12  to  express  the  bending  moment,  the  stress  through- 
out  the  thin  plate  is  given  by  a   ■  6M/H  . 

For  -b  • .  x  *      -c 


6D  r 

a..  -  -p  L m  (A 

H 


cos  mx  +  A,  sin  mx 


m  D 


(35) 


For  -c  <  x  <  0 


-6d 

2    H2 


-a,x 


B1[(af 


b2)  sin  b2x  -  23^  cos  b2x]' 


B2  [(a2  -  b|)  cos  b2x  +  2axb2  sin  bgx]  + 

p       O 

B_  [(a-  -  b.  )  sin  b.x  +  2a  b  cos  b..x  J  + 
B^[(a2  -  b2)  cos  bxx  -  2a1b1  sin  ^x]  +  0(Hv)  & 


V 


For  0     <  x   <    d 


;6D 
H2 


"alX 


C1[(a1   -   b2)    sin  b?x  -   23^  cos  b2x  ]  + 


2  2 

C2  ^ai   "   b2^   C°S  b2X  +  ^l^  Sin  b2X^ 


p  p 

C     [(a      -   b   )    sin  b  x  +  2a   b     cos   t>,xj  + 


ck  L(ai  ■  bi)  cos  bix  -  2aibi  sln  bix3 


+    P(l-H.)f 


(36) 


(37) 
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For  d  <  x  <  (a-b) 

o.  =  &  [  m2  (D_  cos  mx  +  D^  sin  «)  -  ^f-  -  P  (**)  ^]  (38) 

4-   „2        j  mi) 

For  the  case  where  A  <  0,  the  deflection  of  the  plate  in  zones 
1  and  k  (Figure  l)  is  given  by  Equations  8  and  11  respectively.  In 
zones  2  and  3,  the  deflection  may  be  expressed  as  follows: 

a,  x  ax         -a..x  .      . 

^  B  qo/K  +  B  e  *     +  B6e  3     +  e     l   (B?  sin  bf  +  Bq  cos  bf)  (39) 


„     =  VK  +  CeK%  +  c/^  +  e"a3X(C7  sin  bf  +  Cq  cos  bf)  («*» 


a,  x.     a^x    -aox 
3 

Using  the  same  boundary  conditions,  Equation  7,  and  following  the 
procedure  previously  outlined,  the  constants  may  be  evaluated  and  the 
stresses  and  deflections  determined. 

To  obtain  the  complete  solution  to  the  problem,  the  method  of 
solution  is  applied  to  the  various  regions  which  develop  (for  example, 
see  Figure  2) . 

Part  II 
The  solution  to  the  governing  differential  equation  (Equation  l)  is 
again  given  by  Equation  3,  therefore,  for  the  case  where  ^,  *,,  and  A^O, 
and  the  problem  is  again  described  by  Figure  9b,  (ij  >  0)  the  deflection 

of  the  plate  may  be  expressed  as  follows: 

a  x  ™^i " 

„   =  qo/K  +  e  1   (AK  sin  b.x  +  A,  cos  b  x)  +  e    (A  sin  b,,x  + 
la       o         ;>      i     u      * 

Ag  cos  b  x) 


101 


a,x 


-a„x 


w 


=  qo/K  +  e   (Bc  sin  b.x  +  B^  cos  b.x)  +  e    (B„  sin  b„x  + 
lb      o        5      lb      1  7      2 

Bg  cos  b2x)  (42) 


V 


-a2x 


w2  =  4o/K1  +  e    (C  sin  t^x  +  Cg  cos  7^)  +  e    (c  sin  y^x.   + 

Cg  cos  r2x)  (43) 


a,x 


-a.x 


w_  =  qo/K  +  e   (D,.  sin  b.x  +  D^  cos  b.x)  +  e    (D„  sin  b„x  + 
3      o        5      lb      1  7      2 

Dg  cos  b  x) 

To  obtain  the  particular  solution  of  interest,  the  following 
conditions  must  be  satisfied: 
(a] 
(b) 


(c 
(d 
(e 
(f 
(g 
(h 
(i 
(j 
(k 

(1 
(m 

(n 


wla  (— > 


Wla  <— > 


o/K 


Wla  <°> 


"lb  <°> 


Wia  (°)  "  Wlb  <°> 

•J,  (0)  -  w^  (0) 

Wlb  (0)  '  Wla  (0)  =  P/D 

Wlb  (V  =  W2  <V 

wib  (V  B  W2  (Li} 

Wlb  (L1}  =  W2  (V 

Wlb  (V  "  W2  (V 


w,  (L2) 
W2  <L2> 


w3  (L2> 


w2  (L2)  -  *»3   (L2) 


Ml 


w2  (L2)  -W3  (L2) 


(kh) 


(^5) 
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(•)   W3  (oo)  =  %/Ko 
(P)   W^  (»)  =  0 

As  x  approaches  -  °°  ,    both  w.   and  w   approach  oo  >  therefore  if 
conditions  (a)  and  (b)  are  to  be  satisfied,  then  A  and  An  must  equal  zero, 
Similarly,  conditions  (o)  and  (p)  imply  that  D^  and  D^-  must  equal  zero. 

Applying  conditions  (c)  through  (n)  to  Equations  1+1  through  kk,    the 
following  equations  are  obtained: 
A6  -  B6  -  Bg  =  0        •  (k6) 

blA5  +  alA6  "  blB5  •  alB6  "  b2B7  +  alB8  =  °  W) 

2alb]LA5  +  (.J  -  b2)  a6  .  ^^   .  (.2  .  fc2}  ^  +  ^^  . 

(a*  -  b2,)  Bg  -  0  (W) 

(bl  "  3biai>  A5  ■  (al  ■  3»ibi5  A6  '  (bl  ■  3blai>  B5  +  (al  ■  3albl)  B6  + 
(3a^b2  -   b3)   B?  +  (38^  -  a3)   Bq  -  P/D  (1*9) 

alLl^   ,  .  .   .  _      ,  .  ,    "alLl 


e    (B  sin  b^  +  Bg  cos  b1L1)  +  e     (B  sin  b A.     +  BQ  cos  b  L-) 

a2Ll  "a2Ll 

e    (C5  sin  7^  +  Cg  cos  7X\)   -  •     (C^  sin  7^  +  Cg  cos  7^) 


K  -  JL 

'"(i-r1)  '     (50) 

1  o 
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a.L 


e  1  1   [B  (bx  cos  b^  +  a1  sin  b^)  +  B6(a1  cos  b^  -  b±   sin  b^)] 


■alLl 


+  e  [BT(b2  cos  b2Ll  '  31  Sin  b2Ll^   "  B8^b2  Sin  b2Ll  +  al  °OS  b2Ll^ 


a2Ll 


-  e  *      [C   (7X  cos  71L1  +  a2  sin  7^)  +  C6(a2  cos  7^   -  y±     sin  7^)] 


■a2Ll 


e     [  CT(^2  cos  7^  -  &2   sin  7^)  -  Cg^  sin  7^  +  a2  cos  7^)] 


=  0 


(51) 


alLl 


2       ,2, 


B5  [(a2  -  b^)   sin  b^  +  2a1b]>  cos  b^]  +  Bg  [  (a^  -  bp   cos  b^ 


"  2a1b1   sin  b^] 


+  e 


"alLl 


B7[(al   "  bl)   sin  b2Ll   ■  2*\h2  C°S  b2Ll] 


+  Bg  [2axb2  sin  b^  +  (a2   -  b|)   cos  b^] 


-e 


a2Ll 


2  2 


2  2, 


C5[(a2  "  71>   Sin  71L1  +  2Vl   C°S  7lLl]+  °6C(a2  "   71)   C°S   71L1 


2a27x  sin  7^] 


"a2Ll 


C?[  (a2  -   7g)    sin  ?2L1   "   2a272  cos   7^  ] 


2  2 

+  Cg[2a272  sin  7^  +  (&2  "  7g)   cos  72L1^ 


=  0 


(52) 
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alLl 


BJ(a^  -  3a1bj)  sin  b^  -  (b3  -  3^3^)  cos  b^  ]  + 


B6  [(b3  -  3blfl2)  sin  b1L1  +  (a3  -  38^^)  cos  b^  ] 


"alLl 


2  _  „3i 


B7  ^(3a?2  '  b2)  C°S  b2Ll  +  (3alb2  *  31)  Sln  b2Ll  3 + 


3  _    ,.2, 


Bg  [(33^   -  ap   cos   b2L1  +  (b|  -   3ajb2)    sin  b^  ] 


-  e 


a2Ll 


C5^(a2  "   3Vl}   Sin  71L1   "   (71   "   37la2}   C°8  71L1  ^+ 


2x       . 


c6  £(7i  "  3ria2>  sin  7iLi  +  (a2  "  sVi*  cos  7iLi3 


"a2Ll 


2      .3. 


C7l(3a272  '   7P   C°S   72L1   +  (3a272  "  V   SiD  72Ll3  + 


2         3 


3        ,  2. 


Cg[(3a272  -  ap   cos   7^  +  (7£  -   3af,72)    sin  7^] 


=  0 


(53) 


a2L2  'a2L2 

e  '   (C   sin  7XL2  +  Cg  cos  7^)  +  e     (C?  sin  7^  +  CQ  cos  72L2) 


e  8l  2(D  sin  b2L  +  Dg  cos  b^  )  -  %(-^—^) 

1  o 


(5»0 


a  L 

2  2 


e  '   E^7!  cos  7^  +  a2  sin  7^)  +  C6(a2  cos  7^  -  7X  sin  7-^2^~\ 


"a2L2 


+  e      Cc7(^2  COS  72L2  "  S2  SiD  72L2^  "  °8^72  Sin  72L2  +  *2  C°S  72L2^" 


"aiL2 


e  LD7(b2  cos  bgL2  -  a1   sin  b^)    -  Dg(b2  sin  b^  +  &1  cos   b^)] 


(55) 
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a„L. 

.  2  2 


C5  t(a2  '  71)  SiQ  71L2  +  2Vl  C°S  71L2^  + 


C6  t(a2  "  71>  C°S  71L2  "  2Vl  Sin  71L23 


+  e  2  2  [^[(a2  -  72)  Sin  ^ 


2a272  cos  72L2]+  Cg  [(a2  -  72)  cos  7^  +  2a272  sin  7^] 


"alL2 


2   ,2 


D7  ^(ai  '  V  Sln  b2L2  "  2aib2  C°S  b2L2^  + 


.2   .2 


D8  ^(al  '  b2}  C°S  b2L2  +  2alb2  Sin  b2L2^ 


=  0 


(56) 


a  L 
2  2 


C5  C(a|  -  3a272)  sin  7^  -  (7*  -  37^)  cos  7^]  + 


[(7j  "  37xa2)  8in  ^iL2  +  (a2  "  3a27l^  C°S  71L2  M  + 


-artL„ 

2  2 


2        _3> 


C7[(3a272  "   7p   cos  72L2  +  (3a272  -  a|)    sin  7^]  + 


2        _3 


3        ,  2_ 


C8[(3a272  -  ap   cos  7^  +  (7g  -   3a272)   sin  7^] 


"alL2 


D?  [(3a^>2  -  b\)   cos  b2L2  +  (3axb2  -  a\)   sin  b^  ]+ 


(57) 


D8  C(3alb2  "  al)  C°S  b2L2  +  (b2  '  ^lV  Sin  b2L2  3J  *  ° 

The  constants  in  Equations  k6   through  57  were  evaluated  with  the 

use  of  a  7090  computer,  and  the  deflections  of  the  plate  were  determined 

from  Equations  (kl) ,    (k2) ,   (^3)  and  (hh) . 
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The  bending  moment  (M)  of  the  plate  is  given  by  Equation  12,  (AX=0), 
and  the  stress  (per  inch  width)  is  expressed  by  a   ■  6M/H  .  For  the  case 
described  by  Figure  9D>  the  stresses  throughout  the  various  zones  are 
as  follows: 

For  x  <  0 


m   -6De 


ax 


'la    B2 


2   ,2 


A  [(a^  -  bZ)   sin  b  x  +  2a  b  cos  0^]  + 


P    P 
A.- £  (a   -  b.  )  cos  b  x  -  2a  b  sin  b  xj 


For  0  <  x  <  L. 

-6D  r  a'x 


'lb    R2 


1   12    2  i 

e     B_  [(a  -  b. )  sin  b.x  +  2a  b.  cos  b.xj  + 

B/-[  (a   -  b.  )  cos  b..x  -  2a  b.  sin  b  x^  + 


-a,x 


1   BT[(a^  -  b2)  sin  b2x  -  2a^>2   cos  bgx  ]  + 

2    2  -\ 

BQ[  {b1   -   b2)  cos  b2x  +  2a1b2  sin  b^J 


For  L.  <  x  <  L 


2    H2 


a_x 


2    2  -i 

C5  ^32  "  71^  Sin  71X  +  2a27l  C°S  71XJ  + 

C6  ^*2  '   7l^  C°S  71X  "  2a27l  Sin  71X^  + 


-a„x 


2  2 


C7[(a     -  7P)   sin  y  x  -  2a  7     cos  7px]  + 


TL   2   '2 
2/2 


2' 2 


Cfi  [2ap7p  sin  7gx  +  (a2  -  72)  cos  7^] 


(58) 


(59) 


(60) 
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For  x  >  Lr 


-6d 

T     =   — — — 

3    t,2 

J         H 


-a.x 


2    2 
D^[(a1  -  t>2)  sin  bgX  -  23^  cos  bgx  ] 


2    2 
D8[(al  "  V  COS  b2X  +  2aib2  Sln  b2X^ 


(61) 


For  the  condition  when  0  £.  L.  2l    -  &■   (see  Figure  9C)  the  deflection 


of  the  plate  becomes: 
a„x 


-a,x 


W 


=  qo/K  +  e   (A  sin  b.x  +  ks   cos  b  x)  +  e    (A  sin  b  x  + 

AQ  cos  b2x)  (62) 


q 
wA  =  Mo 


a2x 


2 


-/^  +  e  '•  (C  sin  y^   +  Cg  cos  7^)  +  e   "  (C^  sin  ^x  + 
CQ  cos  72x)  (63) 


2b 


a2x 


-a2x 


o/^  +  e   (C  sin  7^  +  C10  cos  7^)  +  e   "  (Cn  sin  y^x.   + 


C12  COS  72X^ 


(6U) 


a,x 


-a.x 


w  =  qo/K  +  e   (D  sin  b  x  +  Dg  cos  b.x)  +  e  1  (D_  sin  b  x  + 
Dg  cos  b  x) 

and  the  boundary  conditions  which  must  be  satisfied  are: 

(a)  Wl  (-00)  =  qo/Ko 

(b)  w'   (-00  )  ■  0 

(c)  w^)  =  w2a(Lx) 

(d)  w;(Ll)  =  w2a(Ll) 

(e)  *lih)   =  w^(Ll) 

(f)  wj^)  -  W-O4) 


(65) 


(66) 
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(g)    w 


2a 


(h)     w' 


(i)  « 

(J)  w 

(k)  w 

(1)  w 

(m)  w 

(n)  w 

(o)  w3 

(p)  «o 


2a 
it 

2a 
in 

2b 

2b 

i 

2b 
ii 

2b 

Ml 

2b 


0)   -  w2b(0) 

o)  =  »'2h(o) 

0)   -  W^CO) 
0)   -  w™(0)   =  P/D 
L2)   =  w3(L2) 
L2)   =  w3(L2) 
L2)  =  w^^) 
L2)   =  w'"(L2) 
oo  )   =  qo/Ko 

oo  )    -    0 


Once  again,    the   deflections  and  stresses  may  be  evaluated  by  applying 
the  boundary   conditions   (Equation  66)    to  Equations  62,   63,   64  and  65  and 
solving  for   the  constants   in  the  resulting  set  of  linear  equations. 

To  solve  the  problem  when  L    <0  (Figure  9a)»    the  deflection  of  the 
plate   is  expressed  as   follows: 

ax  ~aix 

qo/K     +  e        (A,_   sin  b.x  +  A^  cos  b,x)   +  e  (A_  sin  b.x  + 

o  v   5  16  1   '  x   7  2 

Ag  cos  b  x)  (67) 


w 


a,-,x 

_  q^/t,    j.  .  2 


"a2x 


w     =  ^/l^  +  e  c  (C     sin  7.x  +  Cg  cos  7jx)  +  e  (C?  sin  7  x  + 

Cg  cos  72x)  (68) 


w 


.  q. 


a,x 


-a,x 


o/K     +  e       (D-  sin  b.x  +  D^  cos  b,x)  +  e     X  (D„  sin  bnx  + 
3a  o  v5  lb  1'  x7  2 

Dg  cos  b2x)  (69) 
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a„x 


-a,x 


w3b  "  *° 


/K  +  e   (D  sin  b  x  +  D   cos  b  x)  +  e    (D   sin  b  x  + 


D12  C°S  V^ 


(70) 


and  the  following  conditions  are  used  to  evaluate  the  constants: 


(a! 

1     «x   ( 

(b] 

1      wx    ( 

(c] 

I      wx    ( 

(d] 

1 

>      wx    ( 

W 

1  wl ( 

a: 

III  . 

'      1 l 

(g: 

1       W        I 

(h] 

1     . 

1      «2    ( 

(i. 

1      w^    ( 

(j: 

III  , 

1     w2( 

(k 

1  W3a( 

(i. 

1     . 

1      w_    ( 

(*: 

ii    . 

1      w_    ( 
3av 

(n 

>     w3b< 

(o 

>  V 

(P 

>    w3b( 

-  oo  )    i 

=  Vk 

o 

-   00  ) 

=  0 

v- 

w2   (V 

Ll>- 

w'    (Lx) 

V- 

W2   Ch) 

V" 

L2>   = 

W 

L2>" 

%<V 

L2)C 

%<v 

L2)  = 

»'/,         X 

W 

0)     = 

w3b(0) 

0)      = 

W3b<°> 

o)    - 

w''b(0) 

o)    - 

in  .    . 

W3a«» 

oo    )     = 

<Vk 

(71) 


P/D 


) 


0 


For  the  case  where  the  value  of  A  becomes  less  than  zero  for  a 
particular  zone,  the  procedure  to  be  followed  and  the  boundary  conditions 
are  the  same  as  previously  described,  but  the  appropriate  solution  to 
the  governing  differential  equation  must  be  used  for  the  respective 
zones  (see  Equations  3  and  6). 
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APPENDIX  III 
Typical  Program  for  Part  I  (see  Figure  1) 

C  SOL.     TO    NON-LINEAR    PROBLEM 

C  DELI     IS    GREATER     THAN    ZERO 

DIMENSION    A(  14,  15),     INOEXd*) 

01  MENS  I  ON    AMI),    CR(l),     BT  I  1  )  ,     VR  (  1  ) 

(DIMENSION    W(W>,     AN0HC2),     SIGNAI42) 

DIMENSION    XI14),    AA(1),     8A(1),    BN(1) 

COMMON    X 

COMMON    AA,»A,BN,E,GM,HN,GNM,HNM,GN1,HN1,QN1,FN1,F1,GI,H1,Q1, 

1BU,B12,D11,012,QV,QK,V,P,0,TEN 
5001    FORMAT     (6F10.2) 
5020    1*1 

READ    (5,5001)    AMI),    CR(I),    *H,    DM,    BH,    TE 

AL-480. 

EY-4OO0O00. 

RHO-.27100O00E-02 

P-125. 

AMU-0. 15 

BE-0. 000006 

F-1000000. 

WRITE     16,5001)     AMI),    CR(1),    MM,    OH,    BH,    AL.1E 

WRITE     (6,1002) 
3001    FORMAT     (7E17.8) 
1002    FORMAT     (/) 

00    90    K!*l,3 

HCOUNT'KI 

H-MH*2.«(HC0UNT-1.  ) 

0Q=H»( 150.9/1728.) 

QK-Q0/AM1) 

D-EY»H»»3/(12.M l.-AMU»«2) ) 

TEM»-BE»TE«( l.*AMU)/H 

MKK-3 

IF(     TE     .EQ.     20.     .AND.    CR(1)     .EO.     1.5)    MKK-4 

IF     (     TE     .GT.     25.     .ANO.    CR(1)     .EQ.     1.5)     PKK-5 

IF     (     TE     .GT.     25.     .ANO.    CR(l>     .EQ.    2.0)     *KK-4 

00    91    KK-2.MKK 

VCOUNT-KK 

V-352.» (VCUUNT-1.) 

E»«HO»H»V»»2/0 

QV»OQ/( V««2«RH0»H) 

VRU |«V/(SQRT(S0RTU.«AMI)»D/HRHO»H)»»2) )  I  ) 

DELI       -16. •(*.•( I. -CR(  I  >  »»2  )*V«< I ) ••8- ( 8.-36  .«CR < I )»»2» 
1Z7.«CR(I)«*4 >»VR(I >••*♦*.) 

BT( I )»  SORT  (SORT  (AM  I)/(4.«Q)) ) 

Bl-F  __ 

B2-2.*VR(I)*VR( I)»BT(I)»BT(I)«F 

6  3- (BT( I ) •  •*>•( (¥R( I) ••♦1-1. )»F 

B<i-(VR(  I)«»2)«(CR(  1I«*2)«FMBTII)»»6>  

EPS-.10000000E-04 

XX-0.02 

CALL  SITER  ( XX, EPS , Bl , 82 , B3 , B4 » 

AAII )-ABS(XX ) 

BN(I  l-BTII )«SQRT(2.»VR( I )»«2*< A A( I )/BT (  I ) )«»2 
1*(2.»VR(  I  )»CR('l  )»BT(  I)  )/AA(  I  )  ) 

WRITE  (6,3001)  H,V,VR(1 ),BT( 1), AAII), BN(1), DELI 

GNl-BNl I )«»3-}. «8N( 1 )«AA( 1 (••2 


Ill 


H.M1  =  3.«AA<  1  )  »BNU)  »«2-AA(  1)««3 
QN1=2.«AA( l)«BNl 1) 
FNl=AA( 1)«»2-BN( 1)»»2 
GN*=GN1«BN(1 )-AA(l )»HN1 
HNK  =  BN(  1  )«HN1*AA(1) »GN1 
IF  (0EL1-O.)  83,83,82 
83   GO  TC  91 

82   BA<I)=BT(I)«SQRTI2.»VR<I)««2*<AA(I)/8T(I)>»«2 
1-(2.«VR« I  )»CR(I  )»BT<  I) )/AA(  I )) 
F1=AA( 1)«»2-BA( 1)»»2 
G1=BA(  l)««3-3.«6A(  1)«AA(1)«»2 
H1=3.«AA(1)»BA(1 )»«2-AA(l)*»3 
Q1=2.»AA(1)»BA( I) 
GM=Gl»BA(l)-AA( 1)»H1 
HM=BA(  1  )»Hl*AA«  1 )«Gl 

C4=P»AA(l)/l0«(AA(l)»«2«(4.«AA( 1)»»2*BN(1 >«»2*3. »BA ( 1 > • »2 )♦ 
KF1-FN1  )»<4.«AAU)»«2-FNl*Fl)/4.  )) 
C3=C4»(-4.«AA(1)«.2*FN1-F1)/(2.«Q1) 
Cl=C4«<4.»AA(l)»»2*FNl-Fl)/I2.»QNl) 
MMK=1 

20.  .AND.  V  .EQ.  704.)  MMK=2 

20.  .AND.  V  .EC.  704.  .ANO.  CR(1)  .EQ.  1.5  .AND. 
MMK  =  1 


MMK  =  3 
.ANC.  CR(l) 
.ANC.  CR(1) 

MMK  =  5 
.ANC.  H  .EQ. 

MMK  =  3 
.ANC.  CRU) 


.EQ. 
.EQ. 


1.51 

2.0) 


MMK 
MMK 


■2 
=  3 


8.  )  MMK  =  6 


.EQ.  2.0)  MMK=2 


IF  (TE  .EQ. 

IF  (TE  .EQ. 
1  H  .EQ.  12.) 

IF  I  TE  .EQ.  20.  .AND.  V  .EQ.  1056.) 

IF  (  IE  .E(J.  30.  .ANO.  V  .EQ.  1056. 

IF  (  TE  .EQ.  30.  .ANO.  V  .EQ.  1056. 

IF  (  TE  .EQ.  30.  .AND.  V  .EQ.  1408.) 

IF  (  TE  .EQ.  30.  .AND.  V  .EQ.  1408. 

IF  (  TE  .EQ.  40.  .AND.  V  .EQ.  1408.) 

IF  (  TE  .EQ.  40.  .AND.  V  .EQ.  1056. 

00  211  II=MMK,9 

SCCUNT=II 

SS1=-(48.*(SC0UNT-1. )«DH) 

SS3=AL*SSl 

APL  =  0. 

AVL=EXP(AA(1 )«SS1)«(C4»(G1«SIN(BA( 1 ) «SS 1 )-H l«COS (8A i 1 > »SS 1 ) )- 
lC3o(Hl«SIN(BA(l)«SSl)+Gl*CGS(BA(l)«SSl))) 

AVR=EXPI-AA( 1)»SS3)»(C4»«HN1»C0S(BN1 1 ) • SS 3) *GN1 «SIN ( BN( 1 ) • SS3 ) )♦ 
1C1«(HN1»SIN(BN(  l)»SS3)-GNl«C0S(BN<  1)»SS3)  )  ) 

AMR=0. 

Bll=QQ«COS(SCRT(E)»SS3)/(E«»2»D)-( (  TEM* AMR )»COS< SQR T< E ) • SS3  )- 
1AVR»SIN(SQRT(E) «SS3 ) /SQR T < E ) )/E 

B12=CQ»SIN(SQRT(E)«SS3)/(E»»2«D)-< ( T EM* AMR ) «SIN( SORT  IE ) •SS3 )♦ 
1AVR«C0S(SQRT(E)»SS3)/SQRT(E) )/E 

Dll  =  QQ«COS(SCRT (E)»SSl)/(E»»2»0)-( ( TEM*AML )«COS( SQR T(  E )• SSI )- 
1AVL«SIMSQRT(E) »SS 1 ) /SQRTI E ) )/E 

D12=CQ»SIN(SQRT(E)»SS1)/(E««2»D)-I < TEM+AML ) »SIN( SQRT (E ) »SSl )♦ 
1AVL«C0S(SQRT(E)»SS1)/SQRT(E) )/E 

QP=10000. 

DO  412  J=l,12 

DO  412  K=l,13 
412  A(J,K)=C. 

X(13)=SS1 

SZ=SINIBA(1)«X(13) ) 

RZ  =  C0S<BA(1) »X(  13)  ) 

TZ=EXP(AA«1)«X( 13) ) 


APPENDIX  III 


112 


OZ  =  EXf>(-AA(  1  )«Xl  13)  » 
SNZ=SIN(BN( 1  )»X(  13) ) 
RNZ=C0S(BN(1)»X(13)) 
RMZ=COS(SQRT(E)»X( 13)) 
SMZ=SIN(SQRT<E)«X(13)) 
A(ltl)=TZ«(F 1«RZ-Q1«SZ) 
A<1,2)=TZ»(Q1»RZ*FI»SZ) 
A(1,3)=DZ»(FN1«RNZ*QNI«SNZ) 
A( 1,4)=DZ«(FN1«SNZ-QN1«RNZ> 
A«2tl)=TZ»(Gl»SZ-Hl»RZ) 
A(2,2)=-TZ»(HI»SZ*G1«RZ) 
AJ2,3)=DZ»(HN1»RNZ*GN1»SNZ) 
A(2,4)=DZ«(HN1«SNZ-GN1»RNZ) 
A(3,1>=1. 
A(3,3)=l. 
A(4,1)=AA(1) 
A(4,2)=BA(1) 
A(4,3)=-AA( 1 ) 
A(4,4)=BN(1) 
A(5,  1)  =  F1 
A(5,2)=0l 
A(5,3)=FNl 
A(5,4)=-UNl 
A(6,1)=H1 
A(6,2)=G1 
A(6,3)=-HN1 
AU,4)  =  GN1 
DO  413  J  =  3,6 
DO  413  K=li4 
JJ=K»4 
413  A( J, JJ)=-A(J,K) 

X(14)=SS3-12.»VC0UNT 

IF  (SSI  .EQ.  -432.)  X(14)=12. 

SX=SINIBA(1)«X( 14)) 

RX=C0S(BA(1)«X(14)) 

TX=EXP(AA(1) »X< 14) ) 

DX=EXP(-AA(1 )«X( 14)) 

SNX=SIN«BN( 1 )»X(14) > 

RNX=C0S(BN(1)«X(14) )  

RfX=C0S(SQRr(E)«X(14)) 

SMX  =  SIN(SQRT(E)«XU4)) 

A<7,5)=TX«« AA(l)»RX-BA(l)»SX) 

A<7,6)=TX»(AA(1)«SX*BA< 1)»RX) 

A(7,7)=-DX»(AA(l)»RNX*BN(l)«SNX) 

A(7,8)=DX«(BN(1 )«RNX-AA(1)«SNX) 

A(7,10)=-l. 

A(8,5)*TX»(F1»RX-Q1»SX) 

A(8t6)=TX»(Fl»SX*Ql«RX) 

A(8,7)=DX«(0Nl»SNX»FNl«RNX) 

A(8,8)=0X»(FNI«SNX-0N1«RNX) 

AI<J,5)  =  TX»(G1*SX-H1»RX) 

A(<3,6)  =  -TX»(H1«SX*G1»RX) 

A(9,7)  =  DX«(HN1»KNX«-GN1»SNX) 

A(9,8)=0X«(HNI«SNX-GN1«RNX»   

A(10f5)=TX«RX 
A(  10,6)=TX«SX 
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A(  10,7)=UX«RNX 

A(10,8)=DX»SNX 

A(  10,9)=-1. 

A(  10,10)=-X<  14) 

A(11,1)=TZ»RZ 

A(  il,2)=TZ»SZ 

A(  11,3)=DZ»RNZ 

A(ll  ,4)=DZ»SNZ 

A(  ll,ll)=-l. 

A(ll,12l=-X( 13) 

A(12,1)=TZ»( AA( i)»RZ-BA( l)»SZ> 

A(12,2)=TZ»( AA(l)»SZ*BA( 1)»RZ) 

A(  12,^)=-0Z» (AA(L)«RNZ+BNll )  «SNZ  ) 

A< 12.4)=DZ*(BN( I >»RNZ-AA(1)»SNZ ) 

A<  12,12)=-1. 

All,l3)=QV-fc»(Dll«RMZ+Dl2»SMZ) 

A(2,13)=SQRT(E»»3)»(Dll»SMZ-D12»RMZ) 

A(6,13)=P/0 

A(7,13)=S(jKT  (E)«(B12»RMX-Bll»SMXH-QV«X<  14) 

A( B,13)=QV-b«(Bll«RMX+Bl2»SMX) 

A(9,13)=SQRT(E«»3)«(B11«SMX-B12»RMX) 

A(  10.13)=B11«RMX*B12»SMX*QV«X(  14)»»2/2.-QK 

A  I  11  ,13  1=011 •RMZ»012»SMZ+0V»X( l3)»»2/2.-QK 

A(l2,l3!=SQRT(E)»(Dl2»RMZ-Cll«SMZ)*QV»X(13) 

00  213  1=1,12 

DO  213  K=l,13 

213  A( J,K)=QP»A( J,K ) 

CALL  CROUT  < A  ,  1 2  ,  1 , 1 3, DETERM , INDEX ) 
DO  214  J=l,12 
K  =  l 

214  X( J)=A( J,K) 

WRIIE  (6,3001)  (X(J),  J=l,14) 

WRITE  (6,1002) 

N=14 

EPS=0.1 

ISV.=  1 

CALL  NCNLIN  (N, X .EPS , 1 SW ) 

WRITE  (6,1002) 

IF  (X(14)  .LT.  481.  .AND.  X(14)  .GT.  -0.5)  GO  TO  400 

WRITE  (6,3001)  S-Sl,  X(14),  H,  V,  CR(1) 

GO  TC  211 
400  DO  IC1  JJJ=1,41 

XCCUNT=JJJ 

Y=SS1*( XCOUNT-1. )»BH 

SY=SIN(BA(1) *Y) 

RY=C0S(BA(1 ) »Y) 

SNY  =  SIN(BN( 1  )»Y) 

RNY=CUS(BN( 1  )»Y) 

IF  (Y-XU3))  150,  150,  151 
150  W(JJJ)=X(11)+X( 12)»Y+D11»C0S(SQRT( E ) «Y ) ♦D12«S I N( SQR T ( E ) • Y ) ♦ 
lQV»Y«»2/2. 

AMCM< JJJ)=-D»E»(D11»(C0S(SCRT(6)»SS1)-CCS(SQRT(E )«Y  )>♦ 
1D12»( SI N( SORTIE )  «SS1 )-S I N ( SQRT < E ) »Y ) )*AML/E) 

SIGMA( JJJ)=6.»AM0M( JJJ)/(H««2) 

GO  TO  101 
151   IF  (Y-0.)  152,152,153 


Ill* 


152 


153 
154 


155 


101 


211 
91 

90 


W(  J 

IEXP 

AMC 

IX(3 

2X(  1 

SIG 

GO 

IF 

w(  J 

IEXP 

AMC 

1X(7 

2X(5 

SIG 

GO 

W(  J 

1QV« 

AMC 

1B12 

SIG 

CCN 

WRI 

WRI 

WRI 

WRI 

WRI 

COK 

CCN 

CCN 

GO 

ENC 


J  J  )  =  C 
(  AA(  1 
MUJJ 
)»(QN 
)»(F1 
MAI  JJ 
TO  10 
(Y-X( 
JJ)  =  0 
(  AA(  1 
MUJJ 
)«(QN 
)»(F1 
MAI  JJ 
TO  10 
JJ)=X 
Y««2/ 
MUJJ 
♦  (SIN 
MAUJ 
TINUE 
TE     (6 


TE 
TE 
TE 
TE 


16 
(6 
(6 
(6 


TINUE 
TINUE 
TINUE 
TC  50 


K*EXP(- 

>»Y)»(X 

)=-0«(E 

l«SNY*F 

•RY-C1* 

J)=6.«A 

I 

14))  15 

K*EXP(- 

)«Y)«(X 

)=-0«(E 

1»SNY*F 

•RY-Q1* 

J)=6.«A 

1 

(9)+X(l 

2. 

)=-0«E» 

(SORTIE 

J)=6.»A 

,3001) 
,1002) 
,3001) 
,3001) 
,1002) 


20 


AA(1)»Y)»(X(4)»SNY*X(3)«RNY)+ 

(2l»SY+XI 1)»RY) 

XP(-AA( 1 )«Y)  •(X(4)»(FN1»SNY-QN1»RNY)+ 

iMl'RNY)  )*EXP(  AA(  1  )«Y  )  •  ( X ( 2  )  • (  Fl* SY*Q 1 »RY  )  ♦ 

SY) l-TEM) 

M0M(JJJ)/(H«»2) 

4,154,155 

AA(1)«Y)«(X(8  )»SNY+X( 7)«RNY)+ 

(6)*SY+X(5)«RY) 

XP(-AA( 1) »Y) »(X(  8)*(FN1«SNY-QNI*RNY)+ 

Nl»RNY) )*EXP(AA(  1  )>Y ) • ( X ( 6 ) • I F 1» SY*Q 1»R Y  )  ♦ 

SY) )-TEM) 

M0M(JJJ)/(H«»2) 

0 )»Y*B11 »COS( SORTIE )»Y)*B12«S INI  SORT (E)»Y)  + 

(B11«(C0S(SQRT(E)«SS3)-CCS(SQRT(E)«Y))* 
)»SS3)-SIN(SQRT(E)«Y) )*APR/E) 
M0MUJJ)/(H«»2) 

(WIJJJI,  JJJ=1,41) 


(SIGMA! JJJ), 
SSI 


JJJ=1,41) 


24 


25 
26 


SUBROUTINE  SITER  ( XX , EPS  ,Bl , B2, B3, B4 ) 

XNB=XX 

XNl=Bl»XNB»»6+B2«XNB««4*B3«XNB«»2-B4 

XN2  =  6.«B1«XNB«5*4.»B2»XN8»«3*2.«B3«XNB 

XX=XNB-XN1/XN2 

IF  (ABS(XX-XNB)-EPS)  26,26,25 

GO  TO  24 

RETURN 

ENC 
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SUBROUTINE  NONLIN  ( N,X , EPS, ISM ) 
C      N  IS  NUMBER  OF  INDEPENDENT  VARIABLES  ANC  EQUATIONS 
C       X  IS  INITIAL  ESTIMATE  OF  ROOT  AND  MUST  HAVE  DIMENSION  N  IN 
C  CALLING  PROGRAM 

C      EPS  IS  ALLOWED  ABSOLUTE  ERROR 
C      ISW  IS  INTERMEDIATE  OUTPUT  SELECTOR 
C  1  FOR  NO  INTERMEDIATE  OUTPUT 

C  2  FOR  INTERMEDIATE  OUTPUT 

DIMENSION  F( 14,  14)  ,  GI14),  DELTI14),  XI14),  REX ( 14)  ,  SI  14 > , 
1BESTI14),  AI14.15),  INDEX!  14) 

COMMON  X  

REPS  =  .01»EPS 

DO  1  1=1,14 

GII)  =0.0  

DO  1  J*l»14- 

1  F(  I  ,  J)=0.0 

SSAEX  =  I.E6  

SSAET  =0.0 
IC=0 

2  ISC  =  1 _____ 

K=0 

3  CALL  EVAL(F.G) 

C      X  IS  IN  COMMON  

DO  17  I=1,N 
17  GII)  =  -Gil) 

DO  30  I  =  It  N  

30  Sill  = GII) 
DM  =  1.0 

NN=N*1  

QP=1000.  ~      "   ~~~ 

DO  987  1=1, N 

DO  987  J=1,N 

987  AU,j)  =  QP»FlI,JT~ 
DO  988  1=1, N 

988  At I,NN)=QP»GII)  

CALL  CROUT  IA,N,1,NN7DETERM7TNDEXT 

00  100  IJK=1,N 

100   DELT  IUK!  =  AIIJK,1)  

lF(ABS(DETERM)-l.E-207V,4,iO 

4  K=K+1 

GO  TO  (5,6,7),K       

5  XIISC)=1.1»  XI I  SCI "~ 
GO  TO  3 

6  X(ISC)=  .9.XIISC1/1. 1  

GO  TO  3 

7  X(ISC)=  l.ltX(ISC) 

ISC=ISC  ♦  1  

IFIISC-NI  8,8,9 

8  K=0 

GO  TO  4  

CALL  PUNTIL,  BEST,  IC,  ICB,  N,  ITRPT) 

C      ZILCH       .  

10  DO  Zl  I»liN 

21  REXII)  =  ABSIDELTI  D/XII  )) 
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SS*CT  =  0.0 
DO  23  1=1, N 

23  SSAET  =  SSAET  ♦  S(I )«»2 
IF(SSAET  -SSAEX)  24,24,26 

24  DO  27  I  =  1,N 
27  BESTU  )  =  X(  I) 

SSAEX  =  SSAET 
ICB  =  IC  ♦  1 
26  DO  11  1=1, N 

IF!REX! I  )-REPS)  11,11,12 

11  CONTINUE 

DO  25  I  =  1,N 

IF  (S( I  )-EPS)  25,25,12 

25  CONTINUE 
GO  TO  15 

12  IC=IC*1 

IF(  ISW-D14, 13,14 

13  IFUC-100)18,18,19 

14  CALL  INTEK(N,IC,REX,S) 
IFUC-100)18,18,19 

18  DO  20  1=1, N 

20  X(  I  )  =X(  I  UDELT!  II 

IF  (X(14)  .LT.  -24.  .OR.  X(14)  .GT.  500.)  RETURN 
DIMENSION  CYCLE!20,5) 
JC  =  M0D((IC-1),5)+1 
DO  40  I  =  1,N 

40  CYCLE! I ,JC)  =  X(I) 
IF  (IC-5)  42,41,41 

42  K  =  IC 
GO  TC  A3 

41  K  =  5 

43  DO  50  J  =  1,K 

IF  (J-JC)  44,50,44 

44  DO  46  I  =  1,N 

IF  (CYCLEU.JC)  -CYCLE!  I, J)  150  ,46,50 
46  CONTINUE 

ITRPT  =  IC-MOD! ( 5+JC-J ) , 5 ) 

L  =  3 

CALL  PUNTU,  BEST,  IC,  ICB,  N,  ITRPT) 
50  CCMINUE 

GO  TO  2 

15  SSREX  =0.0 
DO  16  1=1, N 
XII)=XI  I)*DELT(  I  ) 

16  SSREX  =SSREX  *REX!I)»«2 
CALL  FINAL(N, SSREX, SSAEX) 
RETURN 

19  L  =  2 

CALL    PUNT!L,    BEST,     IC,     ICB,    N,     ITRPT) 

RETURN 

ENC 
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SUbROUTINE  PUNT (L, BE ST, IC, ICB,N, ITRPT) 

DIMENSION  X(14>,  BEST(1<.) 

COMMON  X 

GO  TC  (  1,2,7) ,L 

1  ICT  =  IC  +  1 
WRITE  (6,3)  ICT 

3  FOPMATUH  31H  SYSTEM  IS  IN  A  SINGULAR  REGION/1H  35H  SINGULARITY  OC 
1CURREO  ON  ITERATION  U/1H  27H  THE  SINGULAR  POINT  FOLLOWS) 

WRITE  16, M   (XII), 1  =  1, N) 
A  FORMAT! 1H  E20.8) 

STOP 

2  WRITE  (6,5)  ICB 

5  FORMATdH  33H  NUMBER  OF  ITERATIONS  EXCEEDS  ICO/  1H  11H  ITERATION 
114,  44H  IS  BEST  ESTIMATE  SO  FAR  AND  IS  GIVEN  BELOW  ) 

WRITE  (6,6)  (BEST( I), 1=1, N) 

6  FORI-AT(  1H  E20.8) 
STCP 

7  WRI TE  (6,10)  ITRPT, IC 

10  F0PMAT(13H1  ITERATIONS  13,  4H  AND  I3.46H  ARE  IDENTICAL  INDICATING 
1  A  CYCLIC  CONDITION.  M3H    THE  BEST  RESULTS  SO  FAR  ARE  GIVEN  OELO 
2W.   ) 

WRITE  (6,6)  (BESTU  )  ,I  =  1,N) 

STCP 

END 
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SUBROUTINE  INTER(N, IC,RfcX,S) 
DIMENSION  X(l<»),  REX(14)i  S(14) 
COMMON  X 
WRITE  (6,1)  IC 

1  FORMAT* IHO  20H  ITERATION  COUNT  IS  13) 
WRITE  (6.2) 

2  FORMAK1H  5X,  5iH  ESTIMATEC  ROOT 
1  ERROR) 

WRITE  (6,3)(X(I),REX(I),Sm,I=l,N) 

3  FORMAT11H  3E20.8) 
REIURN 

ENC 


RELATIVE  ERROR 


ABSOLUTE 


SUBROUTINE  F I NAL ( N, SSREX ,SS AEX ) 

DIMENSION  X( K) 

COMMON  X 

WRITE  (6,1) 

FORMATdHO  20X,  12H  FINAL  ROOT  » 

WRITE  (6,2)  (X(I),I=l,N) 

FORMATI1H    20X,     E20.8I 

WRITE     (6,3)     SSREX 

FORMAT* 1H  10X,  38H  SUM  OF  SQUARES  OF  RELATIVE  ERRORS  IS  E20.9  ) 

WRITE  (6,4)  SSAEX 

FORMATdH  IOX,  38H  SUM  OF  SQUARES  OF  ABSOLUTE  ERRORS  IS  E20.9) 

RETURN 

END 
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CRCUT  REDUCTION 

SUBROUTINE  CROUT  ( A,N,M, NN, DET ERM,  INDEX  ) 

DIMENSION  A! 14,15) ,  INDEX! 141 

DET=1.0 

JZ=N-1 
JA=N*1 

DO  30  1=1, N 

30 

INCEX! I  )  =  I 

00  700  J=1,NN 

DO  800  11  =  1,  N 

SUP=0.0 

I  =  INDEX(  II  ) 

IFIII-J)33,34,34 

33 

IF  « I  I -1 )  9000,9200,9000 

9000 

LLLL=I 1-1 

DO  910C  K=1,LLLL 

IPPP=INOEX!K) 

91C0 

SUN=SUM*AII,K)»AIIPPP,J) 

9200 

A!  I  ,  J)= ( Al  1  ,  Jl-SUMI/A!  I, II  ) 

GO  TO  8C0 

34 

IF(j-l)  8000,8200,8000 

8000 

LLLL=J- 1 

DO  8100  K=1,LLLL 
IPPP=INDEX!K) 

81C0 

SUK=SUM*A!l ,K)«A!IPPP,J) 

8200 

All , J)=A(I,J)-SUM 

800 

CONTINUE 

IF! J-N)41,700,700 

41 

L=INDEX!J) 

KA  =  L 

HIGH=A(L,J) 

KZ  =  0 

DO  35  l=J,JZ 

JC=I*1 

L=INDEX!  JO 

IFJABS  IHIGHI-ABS  ( A! L , J ) ) ) 36, 35, 35 

36 

HIGH=A(L,J) 

KA  =  L 
KZ  =  l 

35 

CONTINUE 

IF(KZ.NE.  0)   DET=-DET 

IFIABS  «HIGH)-l.E-05)3l,31,3200 

31 

WRITE!6,32)  HIGH 

32 

F0RMATI48H0THE  PIVOT  ELEMENT  IS  LESS  THAN  l.E-05  VALUE  IS  E20.8) 

3200 

DO  37  K  =  1,N 

KK  =  K 

IF!  INDEX (K)-KA) 37,38,37 

37 

CONTINUE 

38 

ITEMP=INUEX! JJ 

INDEX! J)=INUEX(KK> 

INCEX(KK)=ITEMP 

700 

CONTINUE 

IF!M>2000, 1000, 2000 

2000 

L  =  N-l 

DO  39  J=JA,NN 
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LL  =  1 

DO  42  K=1,N 

IMABS  (A(K,  J))-0.0)43,42,43 

42  CONTINUE 
IZ=INDEX(N) 

IF(A6S  (A(IZ.N) )-1.0E-02)46,46,44 

44  WRITE(6,45) 

45  FORMAT! S9H0                                  ONLY  SOLUTION  IS  ZERO  VEC 
1T0R)  

GO  TO  10 

46  A(  IZ,J)=5.0000 
IZZ=IN0EXIN-1) 

IF(ABS  <A(IZZ,N))-1.0E-04)47,47,43 

47  A(IZZ,J)=2. 50000 
LL=2 

43  DO  40  IJ=LL,L  

SUM=0.0 

II=N-IJ 

1= INDEX (II) 

LL=I 1*1 

00  9300  K=LLfN  

IP=INDEX(K) 
9300  SUM  =  SUM1*A<  I,K)«A(  IP.J) 

A(I,J)=A(I, J)-SUM1 

40  CONTINUE  "  "   

39  CONTINUE 
1000  DETERM=1.0  

00  900  1=1, N 

K=IN0EX(  I) 


900  DETERM=DETERM»A(K,I) 
DETERM=DETERM«DET 
DO  400  1=1, N 
DO  400  J=JA,NN 
K=INDEX(I  ) 
L=J-N 

400  A( I,L)=A(K,J) 
10  RETURN 
END 
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SUBROUTINE  EVAL«F,G) 

DIMENSION  F<14,14),  G(14),  AA(1),  BA(1),  BN(l),  X(14) 
COPMCN  X 

COMMON  AA,BA,BN,E,GMtHM,GNM,HNM,GNl.HNl,QNl,FNl,Fl,Gl,HltQl, 
iBll,B12,Dll,012.CV,OK.V,P,D,TEM 
00  512  J=l,14 
00  512  K=l,14 

512  F(J.K)=0. 
RMZ=C0S(SQRT(E)»X(13)> 
SMZ=SIN(SQKT(E)»X( 13  » ) 
SZ=SIN(BA«1 ) »X<  13) ) 
RZ=C0S(BA(1)«X( 13)) 
TZ=EXP< AA(1)»X( 13) ) 
DZ=EXP(-AA(1  )«X(  13)  ) 
SNZ  =  SIN(0N(1  )*X(  13)) 
RNZ  =  COS(BN( 1  )«X(13>) 
F(l,l)=rZ«(FJ»KZ-01»SZ) 
F(1,2)=TZ«(01«RZ+F1»SZ) 

F(  1,3)=DZ»(FM«KNZ*QN1«SNZ) 

F(1,4)=DZ»(FN1»SNZ-QN1»RNZ) 

F«2,1)=TZ«(G1»SZ-H1»RZ) 

F(2,2)=-TZ»(H1»SZ*G1«RZ) 

F(2t3)=OZ»«MNl»KNZ*GNl»SNZ) 

F(2,4)=DZ»(HN1»SNZ-GN1«RNZ) 

F 12. 13 )=T2*( X (1 )«(GM»RZ*HM«SZ)-X(2)«<HM»RZ-GM»SZ)>* 
1DZ»<  X(3)»(GNK»RNZ-HNM»SNZ)*X(4)«1HNM»RNZ»GNM«SNZ  )>- 
2E««2»(D11«RMZ+D12»SMZI 

F(3,l)=l. 

F I  3 ,3  )=  1  . 

F<4,1)=AA<1) 

F(4,2)=bA(l) 

F(4,3)=-AA(1) 

FU,4)  =  BN(1) 

F(5,1)=F1 

F<5,2)=01 

F(5,3)=FN1 

FJ5,4)=-CNl 

F(6,1)=H1 

Fl6,2)=Gl 

Flfc,3)=-HN1 

FU.4 )=GN1 

DO  513  J=3,6 

DO  513  K=li4 

JJ=K+4 

513  F( J, JJ)=-F( J,K) 
SX=SIN(UA(1) «X< 14) ) 
RX  =  COS(BAU)«X(  14)  ) 
TX=EXP( AA(1)«X<14)) 
DX=EXP«-AA< 1 )«X(14) ) 
SNX=SIN(BN«  1  )*X(  14)  ) 
RNX=C0S«BN(l)»X(14)) 
SMX=»SIN(SQRT  (E)«X(14)  ) 
RMX=COS(SURT (E) «X(14)) 
F(7,5)=TX»<AAI1)«RX-BA< 1)«SX) 
F(7,6)=TX«<AA(1 )«SX*BA(1)«RX) 
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F(7, 
F(7, 
F<7, 

F(e, 

F(fi, 
FIB, 
F(8, 
F»9, 
F<9, 
F(9, 
F(9, 
F(9, 

1DX»< 

2E»«2 
F(10 
F(10 
F(  10 
F(10 
F(10 
F(ll 
F(  11 
F(ll 
Fill 
FU2 
F{  12 
F(12 
F112 
Fill 
F«  12 
F(13 
F(13 
F(13 
F(  14 
F<  14 
F(14 
Gil) 

1012* 
G(2) 

1SQPT 
G(3) 
G(4) 

1BA(1 
G(5) 

1FN1« 
G(6) 

1HN1» 
G(7) 

1XJ10 
G(b) 

ie»  <b 

G<9) 

ISQRT 

G(  10 

G(  11 

G(  12 

1SQRT 


7)=-UX»< A 

8)=0X«(BN 

10)=-1. 

5)=TXMF1 

6)=TX»(F1 

7)=DX«(3N 

8)=l)X»(FN 

5)=TX*(G1 

6)=-TX»(H 

7)=0X»(HN 

8>=0X«(HN 

14)=TX«(X 

X(7)«(G^M 

•(B11«RMX 

,5)=TX»RX 

,6)=TX»SX 

,7)=DX«RN 

,b)  =DX»SN 

,14)=F(7 , 

,1)=TZ»RZ 

,2)=TZ»SZ 

,3)=OZ«RN 

,4)=DZ«SN 

,1)-TZ»( A 

,2)-=TZ»(  A 

,3>=-DZ»( 

,4)=DZ«(B 

fl3)=F(12 

,12)=-1. 

i9)=l. 

,10)=X(14 

,14)=X<10 

,111=1. 

,12>=X( 13 
,13)=X<12 
=F( 1,1)»X 

SMZJ-OV 
=F(2,1)«X 
(E«»3)»(C 
=  X(  1KX(3 
=  AAU)»X( 
)  .X<  M+AA 

=Fi.xm* 

X(7)*QN1« 
=H1«X(1  )♦ 
X(7)-GNl» 
=F(7,5)«X 
)-SQRT(E) 
=F(8,5)»X 
11«RMX+B1 
=F(9,5)»X 
(E»«3)»(B 
)=F(10,5) 
)  =  F (11,1) 
)=FJ 12,1  ) 
IE>*(012« 


A<  1 )«RNX*BN(i) «SNX) 
<i  )«RNX-AA< 1)»SNX) 

•RX-Q1«SX) 

•SX+Q1»RX) 

l«SNX+FNl«RNX) 

1»SNX-0N1«RNX) 

•SX-H1»RX) 

1«SX*G1«KX) 

l»RNX+GNl*SNX) 

1»SNX-GN1»RNX) 

<5)»(3M»RX*HM»SX )-X I  6  )  • ( HM»RX-GP»SX ) )♦ 

•RNX-HNM»SNX)+X(8)»(HNM«RNX+GNM»SNX) J- 

♦B12»SMX) 


X 
X 
5>»X(5)*F(7,6)»X(6)*F(7,7)«X(7HF(7,8)«X<8) 


Z 

Z 

A( 1)«RZ-BA(1)»SZ ) 

A( 1)»SZ+BA(1)»RZ) 

AA(1)»RNZ*BN( 1 )«SMZ) 

N«1)«RNZ-AA«1)«SNZ ) 

,1)«X(1)*F(12,2)»X(2)+F(12,3)«X(3)*F(12,4)»X<4) 


) 

>-SQRT(E)»(Bll»SMX-B12*RMX)+QV»X( 14) 

) 

)-SQRT(E)«(Cll»SMZ-012»RMZ (♦0V*XI13) 

(1 >+F( 1»2)«X(2)*F( l,3)«X(3)-»F(l,4)«X(4)*E*(Dll«RMZ  + 

(l)*F(2,2)«X(2)*F<2,3)»X(3>*F(2,4)»X(4>- 

11«SMZ-012»RMZ) 

)-X(5)-X(7) 

l)+eA( 1)»X(2)-AA( 1  )»XC3)*BN( 1 ) • X ( 4 )-AA ( 1 ) »X ( 5 )- 

(1  )»X(7)-BN(1).X(8) 

Q1»X«2)*FN1»X(3)-QN1»X(4)-F1»X(5)-Q1«X(6)- 

X(BI 

G1»X(2)-HN1»X(3)*GN1»X<4)-H1.X«5)-G1»X(6)* 

X<a)-P/0 

(5)*F(7,6)«X(6)*F(7,7)»X(7)*F(7,8)»X(8)- 

•(B12»RMX-Bli«SMX)-QV»X« 14) 

(5)*F(8,6)«X(6)+F(8,7)«X(7)*F(8,8).X(8)* 

2»SMX)-QV 

(5)+F(9,6)»X(6)*F<9,7)«X(7)*F<9,8)»X(8)* 

12»RMX-Bll»SMX) 

•X<5)*F(10,6)»X(6)*F« 10,7) «X( 7  )  *F  <  10  ,8  )  »X «  8  )*QK 

•  X(l)*F(ll,2)»X(2)*F(ll,3)«X(3)*F(ll,4)»X(4KQK 

•  X<l)*F(l2,2)«X(2)*FU2,3)«X(3)*F(12f4)»X(4)-Xll2)- 
RMZ-Dll»SMZ)-GV»X< 13) 
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rM*l-xl<mx(10»«xm)*Bll«RMX*Bl2»SMX»CV»XCW)«»2/2. 
G!ir»:;un"xU3,ix!l2».Dll.RMZ*D12.SMZ>QV.X(13).-2/2. 


FCltl3»»GI2) 
F(12,13)=G(1  ) 
F(7,14)=G«B) 
F(8,1<»)=GC9) 

RETURN 
ENC 
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C      PVM  CAL.  OEFL.,  MOMENTS  AM:  STRESSES  IN  FULLY  SUPPORTED  SLAB 

C      THfc  SLAB  IS  SUBJ  TO  MVNG  LOADS  AND  HAS  CIFF  K-VALUES 

C      ALL  CEL  GREATER  THAN  ZERO 

C      WW  AND  SIGMX  SHOULD  BE  DIMENSIONED  ACCORCING  TO  NN 

DIMENSION  A( 12,13),  INDEXI12).  DELI2),  VWI500),  SIGPA(500) 

DIMENSION  AA(3) ,eA(3),BN<3) ,W0( 3),VR(3)  ,BT(3),AK(3)  ,CR<3),  AM0(3) 
5001  FORMAT  (7F10.2) 
5020  READ  (5,5001)  (AMI),  1  =  1,2),  (CR(I),  1*1,2),  ALEN,  6H,  AH 

AML  =  C15 

EY=*C0C0O0. 

RHC=.27lOOOOCE-02 

P=125. 

DO  8C  KKK=1,* 

HCCUNT=KKK 

H  =  6.* (hCUUNT-1.  I  »2. 

QU=H«(150. 9/1728.) 

D=EY«H«»3/( 12.»(l.-AMU»»2)  ) 

QP=.10000000E*03 

IF     (H     .EO.    6.)     UP=.10000C0CE*08 

1001  FORMAT  (*E18.8) 

1002  FORMAT  (/> 

1003  FORMAT  (8E1*.6/8E1*.6) 
100*  FORMAT  (9E13.5) 

DO  91  KK=1,5 

VCCUNT=KK 

V=352.»(VC0UNT-1.I 

DO  92  1=1,2 

BT( I )  =  SORT(SCRTI AM  I )/(*.»0)  )  ) 

WOII  )=P«BT<  I  )/(2.«AKI  I)  )*Q0/AM  I  ) 

AMC(  1  >=P/(*.»BT(  I)  ) 

VR(I )=V/(SORT(SORTI*.«AK(I )»D/< (RH0«H)««2)  )  )  )  

DELI  I  )=16.»(*.«(1.-CR<  I  )»»2)«VR(  I  )  ••8-  (  t>.  -3  6.«CR (  I  )«»2* 
127.»CR ( I )••* )«VR(I )••*♦*.) 

WRITE     (6,100*)     UEL(I),V,    H  

IF     (DEL(I)-O.)     81,81,93 
81       WHITE     (6,1003)    AMI),    AK(2),    CR  (  1 )  ,    CR(2) 

GO    TO    8* 
93      F=.10000000E*12 

B1  =  F 

B2=2.«VR(I)»VR(1)«BT(I)»BT(  I  )  «F 

B3=(BT(I)»»*)«((VR(I)»«*)-1.)«F 

B*=(VK(I)»»2)»(CR(I)«»2)«F«(BT(I)»»6) 

X=.02 

EPS=.1COOOOOOE-06 

CALL  SITER  < X ,E PS ,B1 ,B2 , B3, K* » 

AA( I )=ABS(X) 

BMI)=BT(I)»SQRr(2.»VR(I)»«2*(AA(I)/BT(I))«»2 
1+(2.«VR(  I)»CR( I  )»BT( I  ) )/AA(  I)  ) 

BA(I)=BT(l)»SQRT(2.»VR(l)««2*(AA<I)/BT(IJJ»»2 
1-(2.«VR( I )«CR(I ) «BT(  I ) )/AA(  I  )  ) 
92   CONTINUE 

SIG  =  AMCm*6./(H"2J  _   

WRITE  (6,100*)  AK(1)  ,AK(2),CR(  I  ),CR(2),  ALEN, AH, BH,  AH.X2 

WRITE  (6,1003)  AA(1 ),AA(2),8A( I ),BA(2),HN( 1)  ,BN(2), VR ( 1),VR(2), 
IBT( l),BT(2) ,W0( 1 ),W0(2)f AMO(l), AM0(2) t  SIG 
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100 


AM3) 
BA(3) 
AA<  3) 
BNI3) 
G1  =  BA 
G2  =  BA 
G3  =  BA 
GNl  =  B 
GN2  =  B 
GN3  =  B 
Hl  =  3. 
H2  =  3. 
H3  =  3. 
HN1  =  3 
HN2  =  3 
HN3  =  3 
Ql  =  2. 
Q2  =  2. 
Q3  =  2. 
QN1  =  2 
QN2  =  2 
QN3  =  2 
F1  =  AA 
F2  =  AA 
F3  =  AA 
FN1  =  A 
FN2  =  A 
FN3  =  A 
AL1  =  1 
00  99 
AL1  =  A 
X2  =  -( 
N=  I  FI 
NN=IF 

00  94 
SCCUN 
S  =  ALE 
SM=S* 
WRITE 
IF  IS 
T1  =  EX 
T2  =  EX 

01  =  EX 

02  =  EX 
TM2  =  E 
0H2  =  E 
0M3  =  E 
Rl  =  CG 
S1  =  SI 
RN1=C 
SN1  =  S 
R2=C0 
S2=SI 
RN2  =  C 
SN2  =  S 
RM2=C 


=AK(1) 

=BA(l) 

=AA(l) 

=eN( 1) 

• 1)»»3- 

(2) «»3- 

I3)»»3- 

N(  1)»»3 

NI2)«»3 

NI3>»«3 

•  AAI  l)» 

•AA<2)» 

•AA(3) • 

.•AA(1 » 

.■AA(2) 

.•AAI3) 

•AA(1) • 

•AAI2)» 

•AA(3»« 

.•AM  I  ) 

••AAI2) 

. «AAI3) 

(1  )"2- 

i2)»»2- 

(3)»»2- 

AI1  )«»2 

A<2)»«2 

AI3)«»2 

2. 

KJK=1, 
Ll»2. 
12.+AL1 
XII. +12 
IX(  1.-2 

1  =  1, N 
T=l 

N-AHMS 
AL1 

(6,100 
-0.)  19 
PIAAI I) 
P< AAI2) 
P(-AA(1 
PI-AAI2 
XP( AA(2 
XP(-AA( 
XP(-AA( 
SIBAIl ) 
NIBA(l) 
0S(BN(1 
INIttNIl 
S(BA(? ) 
N(BA(2) 
0S(bN(2 
INIBNI2 
0S(BA(2 


3.« 

3.» 
3.* 

-3. 
-3. 
-3. 

6A( 
BA( 
BAI 

•  BN 

•  BN 

•  BN 
BAI 

eAi 

BAI 

•  BN 

•  BN 

•  BN 
RAI 
BAI 
BAI 
-BN 
-BN 
-BN 


UA(l) 
BAI2) 
BAI3) 
•BNI  1 

•  6NI2 

•  BNI3 
l)»»2 
2)»»2 
3)»»2 

I2)«» 

I3)»« 

1) 

2) 

i) 

(1) 

(2) 

(3) 

1)«»2 

2>««2 

3)»«2 

!!)•• 

I2)»» 

I3)»» 


•  AAI  1)«*2 
•AA12>««2 
•AAI 3I«»2 
).AAI 1) ..2 
)*AA(2) *«2 
).AA(3 ).«2 
-AAI 1).«3 
-AA(2)»»3 
-AAI3)««3 
2-AAIl>»»3 
2-AAI2>»»3 
2-AAI3)»»3 


••ALEN*AL1»/AH» 
.•X2/BH) 


COUNT-1. ) 


4)  S,  SM, 
9,100,100 
•S) 

•  S) 
).S) 
)»S) 
)»SM) 
2)«SM) 
3>»SH) 

•  S) 

•  S) 
)*S) 
)»S) 

•  S) 

•S)  _ 
I  »S1 
l»S  1 

>»SM> 


AL1 


126 


101 


102 


SK 
RN 
SN 
RN 
SN 

ou 

00 

A( 

A( 
A( 
A( 
A( 
A( 
A( 
A( 
Al 
A( 
A( 
A( 
A( 
A( 
A( 
DO 
KJ 
00 
A( 
A( 
A( 
A( 
A( 
A( 
Al 
A( 
A( 
A( 
A( 
A( 
A( 
Al 
A  ( 
Al 
Al 
A  I 
A  I 
A( 
Al 
A  I 
Al 
Al 
Al 
A( 
Al 
Al 
A( 
A  I 
Al 


2  =  SIN 
K2  =  C0 

*2=si 

*3=C0 
C3=SI 
101 
101 
IJtJI 
1,4)  = 
1,6)  = 

2.3)  = 

2.4)  = 

2.5)  = 

2.6)  = 

3.3)  = 

3.4)  = 

3.5)  = 

3.6)  = 
A, 3)  = 

4.4)  = 

4.5)  = 

4.6)  = 
102 

=  KI-4 

102 
KL.KJ 

5.3)  = 

5.4)  = 

5.5)  = 

5.6)  = 

5.7)  = 

5.8)  = 

5.9)  = 
5,10) 
fc,3)  = 

e,4)  = 

6.5)  = 

6.6)  = 

6.7)  = 
fc,8)  = 
6,9)  = 
6,10) 

7.3)  = 

7.4)  = 

7.5)  = 

7.6)  = 

7.7)  = 
7,81  = 
7,9)  = 
7,10) 

8.3)  = 

8.4)  = 

e,5)  = 

8.6)  = 

8.7)  = 


(BA(2  )»SM) 
SIBNI 2)«SPl 

N(BN(2)«SK) 
SIBN(3)«SC) 

N(BN( 3)»SH) 

IJ=1,12 

JI=1,13 

)=0. 

-I. 

-1. 

-BNU 

AA(l) 

-BA(  1 

-AA(1 

ONI 

-FN1 

-01 

-Fl 

-GN1 

HN1 

-Gi 

-HI 


KI=5,6 


KL  =  1, 
)=-A( 
01«SN 
D1«RN 
T1«S1 
T1«R1 
-U2»S 
-D2«K 
-T2«S 
=  -T2» 
D1»IB 
-01»( 
T1«1A 
T1»(A 
D2»(  A 
D2»(A 
-T2»( 
=  T2«< 
Di«<F 
D1»(F 
T1MF 
T1.(F 
D2«<0 
-D2»( 
-T2»< 
=  T2»( 
Ol«(H 
01»(H 
-TIM 
T1«(G 
D2«(G 
-D2»( 


4 

KL.K1) 

1 

1 


N2 
N2 

2 

R2 

Ml 

AAI 

All 

AIL 

A(2 

A(2 

AAI 

BAI 

M» 

Ni» 

l»S 

l»R 

N2» 

FN2 

F2« 

Q2» 

M« 

Nl» 

HI* 

l»S 

N2» 

HN2 


)  •  RN  l_nA  A  LI  >«SN  1 )_ 

1  )»KNL*8N(  1  )«SN1  I 

) <S1  +  BA(  1  )•«!  ) 

)«R1-BA<1)«S1J 

)»SN2-BN(2)«RN2I 

)«RN2*BN<2)«SN2) 

2)»S2*BA(2>«R2I_ 


2)»S2-AA(2)«R2) 
SN1-0N1«RN1) 
KN1+QN1»SNU 


1*QI»R1) 

1-Q1»SI) 

RN2-FN2*SN2) 

•RN2*QN2»SN2> 

S2*Q2«R2) 

S2-F2-R2) 

SM-GN1«RN1> 

RN1+GN1»SN1) 

S1»G1»R1) 

1-H1*R1I 

BN2-HN2-SN2) 

•RN2*GN2«SN2> 
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A(e,  9)=T2»(G2«R2*H2»S2> 

A<8,10)=T2«(H2«R2-G2»S2) 

A(9,7)=DH2*SNM2 

A(9,8)=DK2»RK.M2 

AI9,9)=TM2»SM2 

A(9,10)=TM2»RM2 

A(9,ll»=-DM3»SNM3 

A(9,12)=-DM3»RNM3 

A(10,7)=DM2» <BN(2)«RNM2-AA(2)«SNH2) 

A(  10,8)=-DM2«< AA(2)»RNM2*BN<2 )»SNM2) 

A(  10,9>=TM2*(AA(2)*SM2*BA(2)*RM2) 

A(10,10)=TM2«(AA(2) •RM2-BA(2)»SM2» 

A(10,ll)=0M3«(AA(3)»SNM3-BN(3)«RNM3) 

A(  10,  12)=DM3.(AA(3)«RNM3*BN<3)«SNM3) 

AI  11,7)=DH2«(FN2»SNM2-QN2»RNM2» 

A ( 11 , 8 )=DM2» ( FN2»RNM2*QN2»SNM2) 

A<  11,9)  =  TM2«(F2*SM2*Q2»RM2) 

A(  11,10>=TM2«(F2«RM2-Q2»SM2» 
A(U,11I=0M3«(QN:»RNH3-FN3«SNM3) 

A(ll,12)=-OM3»(UN3«SNM3*FN3»RNM3) 

A(  l2,7)=nM2*(HN2»SNM2-GN2««NM2) 

A(  12.8)=DM2»(HN2»RNM2*GN2«SNM2) 

A(12,9)=-TM2«(H2«SM2*G2«RM2) 

A (12,  10)=TM2«(G2»SM2-H2»RM2> 

A(12,11)=DM3«(GN3»RNM3-HN3»SNM3) 

A(  12,12)=-DM3«(HN3»RNM3*GN3»SNM3) 

AI«,13)=P/D 

A(5,13)=GQ»<  AMD-AM  2)  ) /<  AM  1 )  »AM  2  )  ) 

A<9,13>=QQ»<AM2)-AM31)/JAM2A»AM3)]_ 
DO    411    KI=l, 12 
DO    All    KL=1,13 
All    A(KIfKL)=CP»A(KItKLJ 


CALL    CROUT    ( A,12,l,13,DETERM, INDEX) 

DO    103    J=l,NN 

XCCUNT=J  


V=X2*(XC0UNT-1. )«BH 
IF     (Y-S)     104,104,120 
IC4    RY=COS(BA(l)«Y)      


SY=SIN(BA(l)«Y) 
TY=EXP(AA(1I «Y) 
DY=EXP(-AA( I )»Y) 


RNY=C0SIBN( l)»Y) 
SNY  =  SIMBN<1 UYI 
IF     IY-0.)    105,105,106 


105    W=TY«<AC1,1)«SY*A(2,1)«RY)*QQ/AM1> 

AMOM=-D»TY»( All  ,1>»<F1«SY*G1»RY)*A(2,1)«(F1«RY-Q1«SYI) 
GO    TO    501 


106    H»TY»|AI5(1I •SY*A(6,1)«RY)*0Y»( A(  3, 1)»SNY*A  <  4,  1>«RNY)*QQVAM  1 1 
AMCM=-DMTYMA(5,1>»<F1«SV*QI»RYI*A<6,1)«<F1»RV-Q1«SV)»* 

1DY«(A(3,  1)«(FN1>SNY-QN1«RNY >*A<4,1 )«(QN1«SNY*FN1«RNV))1 

GO    TO    501 

120  IF    (Y-SM)    121,131,131 

121  RY=C0S(BA(2I»Y) 

SY'SINiBAI2)«Y) 

TY  =  EXP(AA(2>  «Y) 

DY  =  EXP(-AA(2  )«Y) 
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RNY=C0S(UN(2 )  «Y  ) 

SNY=SIMBN<2  )»Y) 

W=IY»( A<  9,1)  »SY*A(  10,1  )»KY )+DY» I A<  7,1)*SNY*A(8,1)«RNY)*CC/AM2) 

AMCM=-OI  rY»(A(9,l)»(F2«SY*Q2«RY)*A<  10,  1  ) • < F 2»K Y-Q2«SY) )♦ 
IDY«(A(7,1 >»(FN2«SNY-QN2»RNY)*A( 8 , 1 ) • (QN2«SNY*FN2»RN Y » ) ) 

GO    TO    501 
131    RNY=C0S(HN( 3)«Y) 

SNY=SIN(HN(3)»Y) 

DY=EXP(-AA(3  )»Y) 

W=CY«(A(11,1 )»SNY*A(12,1)«RNY)*CQ/AM3> 

AMCM=-0»OY«( A(ll ,1)«(FN3«SNY-QN3«RNY  )*A(12, 1 )• IQN3» SNY*FN3«RNY )  ) 
501  WW(J)=H 

SIGMAU)  =  AM0P»6./IH»«2) 
103  CONTINUE 

MAX  =  1 

UHAX-WW(l) 

DO  20  J  =  2,NN      . 

IF  <HW(J)-WMAX>  20,20,10 
10   MAX=J 

WMAX=WW(J) , 

20   CCNTINIE 

RWC AX=WMAX/W0(1 ) 

YW  =  X2«-FLGAT(PAX-l|  «BH 

BETYW  =  HTU)*YW 

WRITE  (6,1001)  YW,  BETYW,  W*AX,  RWMAX 

MAX=1  

SIGMX=SIGMA(1) 

00  30  J=2,NN 

IF  (SIGMA(J)-SIGMX)  30,30,40 
40   MAX  =  J 

SIGMX=SIGMA( J) 
30   CCMINLE 

RSIGX=SIGMX/SIG 

YS=X2*FL0ATIKAX-1)*BH 

BETYS=BT(1)»YS_^    _ 


WRITE  (6,1001)  YS,  BETVS,  SIGMX,  RSIGX 
GO  TO  94 

199  IF  (SM-O.)  300,  300,2_0Q_ . 

200  T1  =  EXP(AAU) «S) 
T2=EXP(AA(2)»S> 
D2=EXP(-AA(2)»S)   


T*2=EXP(AA<2 )»SM) 
OM2=EXP(-AA(2)»SHI 
DM3=EXP(-AAI 3)«SML 
Rl=COS(BA(l)«S) 
Sl=SIN(OA(l ) «S) 

R2=C0S!BA(2)»SI 

S2=SIN(BA(2)«S> 
RN2=C0S(BN(2)«S) 
SN2  =  SIN(BN(2  )  «S )  _ 
R«2=C0S(BA(2 )»SM) 
SM2  =  SIMBA(2)»SN> 
RMf2  =  CCS(BN<2)«Sf_l 
SNK2=SIN(BN(2)»'JM) 
RNM3=C0S(BN( 3)«SM) 
SNF3=SIN(BN{ 3)»SK) 
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00  201  IJ»1, 12 
00  201  J  I  =1 , 1 3 
201  A(IJ,JI)-6. 
A(l,4)=l. 

A<1.6)=1.  

A(2,3)=BA(2> 
A(2,4)=AA(2> 

A(2,5»=BN(2) 

A(2.6)=-AA(2 ) 
A(3,3)=Q2 

A(3tA)=F2  

AO,5)  =  -QN2 
A(3,6)=FN2 

A(4,3)=G2  

A(4.4)=H2 
A(4,5)=GN2 

A(4,6)=-HN2  

DU  202  KI=3,6 

KJ=KI»4 

00  202  KL=1,4_ 


202  A(KL,KJ)=-A(KL»KI) 
AI5.1)»U*S1 
A(5,2)=T1»RI 


A(5,3)=-T2»S2 
A(5,4)=-T2«R2 
A(5,5)*-D2»SN2_ 


A(5,6)=-D2»RN2 

A(t,l)=ri.(AA(l).Sl+BA(l).Rl> 
A ( 6  ,  2  )^J 1 • ( A  A(l >«R1-BA(1)«S1) 


A(6,3)=-T2»( AA( 2 )*S2*BA ( 2 I «R2 ) 

A(6,A)>T2«(BA(2)«S2-AA(2)*R2) 

A(fc,5)=D2»(AA«2»«SN2-BNJ2)«RN2JL 

A(6,6)=D2»(4A(2 )  »RN2*BN < 2 > • SN2 ) 

A<7,U  =  Tl»(Fl»Sl*Ql»Rl» 

A(7,2)  =  T1«(F1 »R I =0 1 • S I ) 


AU,3)=-T2«(F2»S2*Q2«R2) 

A(7,4)=-T2»<F2»R2-Q2»S2) 

A ( 7 , 5 ) ±^02 • (FN2»SN2-QN2»RN2> 

A(7,6)»-D2»l FN2»RN2+QN2«SN2 ) 

A(e,l)=-TI»(H1»S1*G1»R1) 

A«8,2)=-TI»(HI»R1-G1»SIJ 


A(8,3)=T2»(H2»S2*G2»R2) 
A(e,<»l»T2«<H2»R2-G2»S2) 
A  (  e  ,  5  )  =  -02»  ( HN2jLSN2r.GN2 •*HZ±_ 

A( H ,6 ) =-02»< HN2«RN2*GN2«SN2  I 
A(9,7)=TM2»SP2 

A<9,8)*TM2»RH2 

A(S,9)*0M2«SNH2 

A(9,10)-DM2«RNM2 

A ( 9 , 1 1 )*-QH3»SNH3 


A(9,12)*-0M3«RNM3 

A I 10t7)  =  TH2*(AA(2)>SM2*8A(2)«RM2) 

A  (  10  , 8  )  rJ  M2j»±AA  (11  •RM2-BA«2>«SM2) 


A(10,9)=-DM2«(AA(2)«SNM2-BN(2)»RNM2) 

At  10, 10) =-0M2«( AA( 2  )»RNM2*BN(2>»SNM2» 
A  i 10 , 1 1 1 =PH3  «  <  AA ( 3 ) >SNM3-BN ( 3  >  »RNH3  » 
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A(  10,12)=DI*n«(AA(3)»RNM3*BN(3)«SNM3) 

A{  11,7>=TM2»(F2«SM2*Q2«RP2) 

A(11,8)*TM2»(F2»RM2-Q2»SM2) 

A(  11,9)=0H2» (FN2»SNM2-QN2»RNM2> 

A(  11,10)=DM2«(FN2»RNM2*QN2«SNM2) 

A(ll,ll)=-DM3»(FN3«SNM3-QN3«RNM3) 

A(ll,12)=-0M3»(FN3«RNH3*QN3»SNM3) 

A(  12,7)=-TM2«(H2»SM2*G2«RM2) 

A  I  12,8)=-TM2«(H2»RM2-G2»SM2 ) 

A(  12,9)=DM2»(HN2«SNM2-GN2»RNM2) 

A(12,10)=DM2«(H*2»RNM2*GN2»SNM2 ) 

A(  12,  11)=-0M3»(HN3»SNM3-GN3»RNM3) 

A(  12,12)=-DM3»(HN3»RNM3*GN3«SNM3) 

A(4,13)=P/D 

A(5,13)=0Q»(  AK(1)-AK(21 I ) /< AK ( l) »AK ( 2  )  > 

A(<3,13)=QQ»(  AK(2)-AK(3)  )  /  (  AM  2  )  »AK  (  3)  ) 

DO    422    K  1=1 ,12  

DO    422    KL=l, 13 
422     AIKI.KL  )=QP*A(KI  ,KL) 

C*LL    C«OUT    ( A,12, I, l3,DETERMt INOEXl 

DO    203    J=l,NN 

XCCUNT-J 

Y-X2-MXCCUNT-1.  J«BH 

IF     (Y-S)     204,204,210 
204    RY=C0S(FjA(1)»YI 

SY=SIN(BA(1)»Y} 

TY=EXP(AA(1) »V) 

W=TY»(A(1,1)  •  SY*A(2,1)«RYMQQ/AK(1) 

AMCM=-0»TY«( A( 1  , 1)  •  (FJ  »SY*Q  l»RY  )  *A  (_2 !.,  I  )  •  (JF  1j»RYtQ1«SY J I 

GO    TO    502 

210  IF    (Y-SMI    211,214,214 

211  RY=CCS(BA(2)  »Y) . 

SY=S!N(BA(2)«Y) 

TY=EXP( AA(2)«Y) 

DY=EXP(-AA(2)«Y|      

RNY=C0S(BN(2)»Y> 
SNY=SIN(BN(2  »»Y) 


I F    (  Y-O.l    21  2, 212, 2 1 3_ 


212  W=TY»(A(3,l)«SY*A(4,l)»RY)*0Y»( A(  5  , 1  )»SNY*A  (  6, 1  >»RNY  I+QQ/AM2  ) 
AMCM=-D»(TY«(A(4,1)»(F2»RY-Q2«SY)*A(3,1)«(F2«SY*Q2«RY)I* 

1 D  Y  •  I  _AJ  6  ,  D>(qN2«SNY  *FN2«RN  Y  )  *A  I  5  .  1 )  •  (  FN2«  SNY-QN2«RNY  ) )  > 

GO    TO"  5  02 

213  W*TY*(A(7,1 ) •  SY*A(8,1)»RV)-»DY«( A (9, 1 ) »SNY*A( 10, I )«RNY) ♦QQ/AM 2> 
ACCM»-D»J  TY»  (  A(  8 ,1 )  •  j  F2  «RY-Q2«S  Y  )♦  A  (  7  ,  1  )«  (  F2»SY*Q2*RY1I» 

10Y«(A(l6,l)»(QN2«SNY*FN2«RNY)*A(9, 1 ) •( Fh2» SNY-QN2«RKY  )  I ) 
GO    TO    502 

214  RNY»C0S(BN(3)»Y1 


SNY  =  SIMBN(3>»Y) 
DY=EXP(-AA(3)«Y) 
W  =  C Y« (Ail  ill) •SNY+A ( 12jJ  )»RNY)*Q0/AK(3) 


AMCM— D«DY«( All 1,1) • (FN3«SN Y-QN3»RNY ) *A (12,1 )•  (QN3«  SNY*FN3»RNYH 
502    WW(J)=W 

S  IGPA  (  J)=AM0H»6./(H»»2) . 

203  CCKTINUE 

MAX  =  1 

WHAX=»WW(  1) 


131 


12 


22 


*2 


32 


30C 


301 


DO  22  J 
IF  (WW( 
MAX  =  J 

WMAX-WW 
COMINU 

RWf4X=W 

YW*X2*F 

BETYW=B 

WRITE     ( 

MAX*1 

SIGMX-S 

00    32    J 

IF    (SIG 

MAX  =  J 

SIGMX^S 

COMINU 

RSIGX*S 

YS»X2*F 

BETYS=B 

WRITE    ( 

GO    TO    9 

Tl=EXPt 

T2-EXPI 

02*EXP( 

TM2-EXP 

0M2=EXP 

TM3=EXP 

DM3»EXP 

Rl-CCSI 

S1»SIN( 

R2-C0SI 

S2=SIN( 

RN2-C0S 

SN2=SIN 

RM2=C0S 

SH2  =  SIN 

RNK2*C0 

SNP2=SI 

RM3=C0S 

SM3-SIK 

RNf3»C0 

SNK3»SI 

DO    301 

00    301 

A( IJ.Jl 

All,  8)* 

A(l,10) 

*<  2,7>  = 

A<2,8)= 

AC2.9.* 

A (2, 10) 

A(3,7) 

A(3,8) 

A  I  3 ,  S  I  > 

A(3*10) 

AU,7>» 


JJ-WMAX)     22*22.12 


(J) 

E  

MAX/WOU) 

LOATI P AX-1)»BH 

T(1)«YW 

6,1001)  YW,  BETYW,  WMAX,  RWHAX 


IGHAt 1) 

*2,NN 


MA(J)-SIGPX)  32,32.42 


IGHAIJ) 
E 

IGMX/SIG     

LOATI KAX-1)»BH 

TdMYS 

6.10011  YSi  BETYS, 

AA(1) .S) 

AAi2)»s» ;_ 

-AA(2)«S) 
(AA(2)*SH) 

I-AAI2l«SMi_ 

( AA(3)»SM) 

l-AAOI*5H) 

BA(l)  «S» ___ 


S1GMX,  RSIGX 


»SM» 


BA(1)*S) 
BA(2)«S) 

HA(2) «S) 

(BN12 )*S t 

(BN(2)«S) 

(BA(2)«SM) 

(BAI2 )«SM) 

S(BN(2)«SK) 

N(BN( 2)»SM) 

(BA(3)«SM> 

(BA( 3I.SM) 

S(BN(3)»SM)_ 

N(BN(  3)< 

IJ=1, 12 

J I  - 1 1 1 3_ 

J»0. 

1. 

•I. 

BAI3J 
AAI3) 
BNC3)  __ 
*-AA(3) 
03 

F3 

-UN3 
*FN3 
G3 


132 


AUt8!*H3 

A(4,9|=GN3 
A(<i,10)=-HN3 
DO    302    KI=9,10 
KJ*KI*2 
DO    302    KL=1,A 
302    A(KL,KJ)=-A(KL,Kl> 
A(5,1)=T1«S1 
A(5,2)=T1»R1 
A(5,3)=-T2«S2 
A(5,4)=-T2»R2 
A(5,5)=-D2»SN2 
A(5,6)=-D2»RN2 

A(6,1)=T1»(AA(1)»S1*BA(1)«R1)_ 
A<fc,2)=Tl«(AA<l)»Rl-BA(l>»Sl) 
A(6,3)=-T2»(AA(2)»S2*BA«2»«R2) 
A(6,4)*-T2«( AA(2»»R2-BA(2)«S2) 
AU,5)=-D2«(BN(2J»RN2-AA(2)«SN2) 
A(6,6)=D2«(AA(2 ) «RN2 + BN < 2 ) • SN2 ) 

A(7,1)=T1«(QI»K1*F1«SII 

A(7,2)  =  U»(F1»RI-Q1«S1) 

A(7,3)=-T2«(C2«R2*F2»S2> 

A(7,<.)=»-r2«(F2«R2-Q2»S2> 

A(7,5)=-D2»(FN2«SN2-QN2«RN2) 

A(7,6)=-D2»(CN2»SN2*FN2«RN2) 

A(8,l  )=-U»(Hl«Sl*Gl»Rl) . 

A(8,2)*-TI«(H1»R1-G1»SI) 
A(8,3)=T2«(H2«S2*G2»R2) 
A(8,4)=T2«(H2»R2-G2»S2) 
A(8,5)=-D2»(HN2»SN2-GN2»RN2  1 
A(8,6)=-D2»(HN2»RN2*GN2«SN2) 

A(9,3)=TM2»SK2 

A(9»4)=TH2*RM2 
A(9,5)=DM2»SNM2 

A(9,6)=DM2«RNM2      

A(9,7)=-TM3»SM3 
A(9,8)=-TM3»RM3 

A(9,9)=-DM3*SNM3 . 

A(9,10)=-DM3«RNM3 

A(10,3)=TP2«(AA(2)«SM2+BA«2)«RM2) 
A  ( 1 0 , 4  )  =*  TM2»  (  AA  «  2  )  »R  M2-BA  (  2JJLSM2  JL 


A(10,5)=DM2»(BN<2)«RNM2-AA(2)«SNM2) 
A  (  10 i 6 )  =  -DM2 • ( A A ( 2 ) • RNM2 ♦ BN ( 2 ) • SNH2 ) 

A(10,7)=-TM3«(AA(3)»SM3*BA(  3)«RM3) 

A(  10,8)=-TM3»(AA(3)»RM3-BA( 3»»SM3) 
A(  10,9)=-DM3 »(bN(3 )»RNM3-AA(3  LSNM3) 
A(  10,10)=OH3»(AA(3)«RNH3+BM(3)«SNH3) 
A(11,3)=TM2»(Q2»RM2*F2»SP2) 
A(11,*)=TM2«(F2«RM2-Q2»SM2) 

A( ll,5)*OM2»(FN2«SNM2-QN2«RNM2l 

A ( 11  ,6)  =  DM2« (QN2«SNM2*FN2»RNM21 
A (11, 7)=-TM3«(Q3«RM3*F3»SM3) 
A(11,8)«-TM3«(F3«RM3-Q3»SM3) 


A(  11,9)=-0M3»(FN3»SNM3-QN3«RNM3I 
A (11 ,10)=-DM3»(QN3»SNM3*FN3«RNM3) 
A(12,3»=-rM2«(H2«SM2*G2«RM2J 
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433 


304 


305 
306 


307 


A( 12,4)=TM2«(G2»SM2-H2«RH2) 

A< 12,5)=-DM2»tGN2»RNM2-HN2»SNM2  ) 

A( 12,6)=UM2»(HN2»RNM2*GN2»SNM2) 

A  (12. 7)=TM3»(H3»SM3»G3»RP3) 

A  I 12,8>=-TM3»IG3»SM3-P3»RM3  > 

A( 12,9)=0M3« (GN3«RNM3-HN3«SNM3) 

A  I 12,10)=-DM3»CHN3»RNM3*GN3»SNM3) 

AI4,13)=P/D 

A(5,13)=QQ»(AK(1)-AM2) ) / ( AK ( 1 ) «AK { 2  )  ) 

AI9,13)=UQ«( AK<2I-AK(3) )  /  <  AK  <  2  )  •  AK  <  3  )  > 

00  433  K  I  =  I  ,  12 

00  433  KL=1,13 

A(KI,KL)=0P»A(KI,KL) 

CALL  CrtOUT  < A , 1 2, 1 , 1 3,0ETERM, INOEX ) 

00    303    J=l,NN 

XCCUNT=J 

Y=X2*»XC0UNT-1. )«BH        

IF     (Y-S)     304,    304,     305 
.RY  =  C0S(BA(1 )  »Y) 

SY=SIN(BA(1)«Y) . 

TY=EXP( AA(1) »Y) 

W=TY»IA(l,l  I  •SY*A(2,1)»RY)*QQ/AK<11 

AMCM=-0»TY»(  A(l,l)»lFl»SY*qi>RY)*A(2.1)«f F1«RV-QI«SY)> 

GO    TO    503 

IF     (Y-SK)     306,306,307 

RY=C0S(BA(2)«Y) 

SY=SIN(BA(2)«Y) 

TV=EXP<AA(2) «Y) 

DY=EXP(-AA(2  )«Y^ __ 

RNY=C0S(BN(2  )»Y) 

SNY=SIM HN(2 )  «Y  ) 

W=TY«(AI3,II»SY*A<4,  1»  »RY  J  *0Y»  (  A(  5  f  1  J«SNY*A  <  6  ,  1  JLfRN  Y  )*0Q/AM2  ) 

AMCM=-0»<TY»IA(<,,1)«(F2»RY-02»SY)*A(3,1)»(F2»SY*Q2»RY)>* 
1DY«(    A(6,l»»(QN2»SNY*FN2»RNY)»A(5,l)«(FN2»SNY-QN2»RNY))) 

GO    TO    503  . _ 

RY=CCS(BA(3) »YI 

SY=SIN(BA(3) «Y) 

TY=EXP(AA(3) «Y) 

DY=EXP(-AA( J)»YJ 

SNY=SIN(BN(3)«Y) 

RNY=COS(BN«i  )»Y) 

IF     (Y-0.)    308,309,309 


308 

309 
503 
303 

11 


W=TY«(A(7,1)«SY*A(8,1)«RYJ*0Y»(A(9,1)»SNY*A(10.1  )»R  NY  l+QQ/AK  I  3  ) 
AHCM  =  -0»(TY«IA(tt,l)»IF3»RY-Q3«SY)*A(7,H»(F3«SY»Q3»HV))» 


IDY«(A(10,1>»<QN3»SNY*FN3«RNY)»A(9,1)»<FN3«SNY-QN3«RNYJ S) 
GO    TO    503 

W  =  DY»<  A(  11,1  )«SNY*AC  12tll*RNY)*QQ/0KO> 

AMCM=-0»DY«(  A(  11,1  )MFN3»SN Y-QN3-RNY )♦* < 12,1) •  (  QN3»  SNY*FN3»RNY)  ) 
WW I J)=W 

S  I  GPA  (  J  J  =  AM0M*6^ /  (  H»»2  J 

CONTINUE 
MAX  =  1 

WMAX=»WW(l» [ 

DO  21  J=2,NN 

IF  (WW(J)-WMAX)  21,21,11 

MAX  =  J 
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21 


41 


31 


94 

99 

91 
84 
80 
95 


WMAX=WW ( J) 

CONTINUE 

KWMX^HAX/rtCIl  ) 

YW=X2«TL0AT(MAX-1)»BH 

HETYW=BT(l)»YW 

WRITE  (6,1001)  YW,  BETYW,  WWAX,  RWMAX 

MAX=1 

SIGMX=SIGMA( 1) 

00  31  J=2,NN 

IF  (SIGMA(J)-SIGHX)  31,31,41 

MAX  =  J 

SIGMX=SIGMA( J) 

CCNTINUE 

RSIGX=SIGMX/SIG 

YS=X2+FL0AICHAX-1)»BH 

BETYS=BT(l)«YS 

WRITE  (6,1001)  YS,  BETYS,  SIGMX,  RSIGX 

CONTINUE 

WRITE  (6,1002) 

CCNTINUE       

WRITE  (6,1002) 
CCNTINUE 

WRITE  (6,100  2) ; \ 

CCNTINUE 
GO  TO  5020 
END 


24 


25 
26 


SUBROUTINE  SITER  ( X , EPS  ,  Bl . B2 , 83, 64 ) 

XNB  =  X 

XN1=B1»XNB»«6*B2«XNB«»4*B3»XNB«»2-B4 

XN2=6.«B1»XNB»»5*4.«B2«XNB««3*2.«B3»XNB 

X=XNB-XN1/XN2 

If-     (ABS(X-XNB)-EPS)     26,26,25 

GO    TO    24  

RETURN 
ENC 
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CRCUT  REDUCTION 

SUBROUTINE  CROUTIA,N,M,NN,CtTERM, INDEX) 

DIMENSION  A(N,NN>, INDEX(N) 

DET=1.0 

JZ=N-1 

JA=N*1 

DO  30  I=l,N 

30 

INCEXI I)=I 

DO  700  J=l,NN 

DO  800  11=1, N 

SU*=C.C 

I  =  INCEX(  II) 

I  F  (  I  I-J)33, 34,34 

33 

IF(II-l)  9000,9200,9000 

900C 

LLLL=II-1 

DO  9100  K=l,LLLL 

IPPP=INDEX(K) 

9100 

SUP=SUP+A(I,K)«A(IPPP,J) 

9200 

A(  I,J)=IA(I,J)-SUM)/A( 1,11) 

GO  TO  800 

34 

IF(J-l)  8000,8200,8000 

8000 

LLLL=J-1 

DO  8100  K=1,LLLL 

IPPP=IN0EX(K» 

eico 

SUf^SUP  +  AU  , K)»AIIPPP, J) 

82C0 

AU,J)  =  AU,J)-SUM 

80C 

CONTINUE 

IF(J-N)41,700,700 
41  L=INCEX(J) 
KA  =  L 

HIGH=AIL,J)  ...  

KZ  =  0 

DO  35  1=J,JZ 

JC»I*1  _____ 

L=INDEXI  JO 

IFIABS  (HIGH)-AbS  ( A ( L  ,  J  )  ) ) 36, 35, 35 

36  HIGH=A(L,J) 

KA  =  L 
KZ=1 

35  CONTINUE 

IFIKZ.NE.  0)   DET=-DET 

IFIABS  IHIGH)-1.E-05)31,31,3200 

31  WKITE(6,32)  HIGH  

32  F0RMATI48H0THE  PIVOT  ELEMENT  IS  LESS  THAN  l.E-05  VALUE  IS  E20.8) 
32C0  DO  37  K=l,N 

KK  =  K  _ 

IF11NDEX(K)-KA)37,38,37 

37  CONTINUE 

38  ITEMP=INDEXI J)     .  

INCEXI J)=lNOEXIKK) 
INCEX(KK)=ITEMP 

700  CONTINUE  _  _ 

IF(M)2000, 1000, 2000 
2000  L=N-1 

DO  39  J  =  JA,NNJ 
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LL  =  1 

00    42    K=1,N 

IFUBS    (A(Kf  J))-0.0)43,42,43 

42  CONTINUE 
IZ=INDEX(N) 

IF(AOS     IA(IZ,N)  )-1.0E-02)46,46,44 

44  WKITE (6t45) 

45  F0RMATI59H0 
1T0R) 

GO    TO    10 

46  A(  IZ,J)=5.0000 
IZZ=IN0EX(N-1) 

IF(A8S     (A(IZZfN)  )-1.0E-04)4  7,4  7,43 

47  A<  !ZZiJ)  =  2.bOOOu 
11  =  2 

43  DO    40     IJ=LL,L 
SUH  =  0.0 

I I=N-I J 
I=INDEXI II) 

LL=  I  I ♦  I  

00  9300  K=LL,N 
IP=INDEX(K) 
9300  SUPl>SUMl+A(ItK)«A(lP(J) 
A(  I,J)=AU,J)-SUM1 
CONTINUE 


ONLY  SOLUTION  IS  ZERO  VEC 


40 

39 

1000 


900 


4C0 
1C 


CCMINUE 

DETERM=l.O 

00    9C0     I=ltN 

K=INOEX(  I) 

DETERM=DETERM»A(K,I ) 

DETERM=OETERH»OET 

00  400  I=ltN 

00  400  J=JA,MN 

K=INDEX( I ) 

L=J-N 

A(I  ,L)=A(K, J) 

RETURN 

ENC 
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DISTANCE   (Feet) 
4 1,2 Ifi 2a 24 20 42 # kp. 


0.0 


H  ■  8  Inches 

H  «=  10  Inches 

H  =  12  Inches 


100 


K  =  200 


K  =  100 


200 


Figure  IV-1       Deflection  at  Ends  of  Slab  versus  Position  of 
Moving  Load  (  C  =1.5,   AT«=30  Deg.,   v=0  mph.) 
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Figure  IV-2   Deflection  at  Ends  of  Slab  versus  Position  of 
Moving  Load  (  C  =1.5,  AT=30  Deg. ,  v=40  mph.) 
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Figure  IV-3   Deflection  at  Ends  of  Slab  versus  Position*  of 
Moving  Load  (  C  =1.5,  AT=30  Deg. ,  v=30  mph.) 
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